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Background 


The formation process of the planets in our Solar System has been the focus of scientific inquiry 
for centuries. It is generally accepted that all of our planets formed from the material in the 
primitive solar nebula. However, the giant planets do not echo the solar chemical composition. In 
recent years researchers have made great advances in reconciling theories with observational data 
obtained from ground based and fly-by detectors. The recent discovery of extra-solar planets 
(Mayor and Queloz 1995, Marcy and Butler 1996) provides us with a larger and more variable 
sample with which we can test our theories as well as challenge our comprehension of established 
models for planetary formation. 

The best candidate model for planet formation is the “core instability” model which proposes 
initial accretion of solid matter until the core is large enough to capture massive quantities of gas 
from the solar nebula. Early work by Mizuno et al. (1978) and Mizuno (1980) consisted of a 
series of equilibrium models that were constructed with solid cores and gaseous envelopes. 
Bodenheimer and Pollack (1986) were the first to construct models from evolutionary 
calculations. Their models assumed that the solid body accretion rate was time-invariant for a 
given evolutionary sequence and was different for different sequences. From these early studies it 
appears that the core instability is capable of explaining the bulk mass properties of the giant 
planets. Pollack et al. (1986) and Podolak et al. (1988) used the model envelopes of 
Bodenheimer and Pollack (1986) to study the ability of the accreted planetesimals to pass through 
the envelopes of the giant planets and reach the core. They found that a significant fraction of the 
planetesimals would have dissolved in the planetary envelope thus enriching them with high-Z 
elements. Thereby, the core instability model is capable of explaining the nonsolar atmospheric 
compositions that are observed. 

Our approach to improving the core instability model was to construct evolutionary models of the 
forming giant planets that allow for the interactions of planetesimals with the envelopes of the 
giant planets, and to calculate the rate of planetesimal accretion rather than to prescribe the rate. 
We selected the accretion model of Lissauer (1987) because it offered a promising means of 
solving the timescale problem and it represented an opposite extreme from the prior assumption 
of a constant planetesimal accretion rate. This provided an opportunity to examine whether there 
is a qualitative difference in the results. 

All simulations are characterized by three major phases. During the first phase, the planetesimal 
accretion rate, which dominates that of gas, rapidly increases owing to runaway accretion, then 
decreases as the planet's feeding zone is depleted. During the second phase, both solid and gas 
accretion rates are small and nearly independent of time. The third phase, marked by runaway gas 
accretion, starts when the solid and gas masses are about equal. The overall evolutionary time 
scale is generally determined by the length of the second phase. 
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We judge the applicability of a given simulation to planets in our Solar System using two basic 
yardsticks. One yardstick is provided by the time required to reach the runaway gas accretion 
phase. This time interval should be less than the lifetime of the gas component of the solar 
nebula, t sn for successful models of Jupiter and Saturn, and greater than tsn for successful models 
of Uranus and Neptune. Observations of accretion disks around young stars suggest that ^ ^ 

1 0 7 years, based on observations of the dust component. The lifetime of the gas component is 
unknown, but it may be somewhat longer (Strom et al. 1993; but see also Zuckerman et al. 1995). 
A second yardstick is provided by the amount of high-Z mass accreted, Mz In the case of Jupiter 
and Saturn, Mz at the end of a successful simulation should be comparable to, but somewhat 
smaller than, the current high-Z masses of these planets, since additional accretion of 
planetesimals occurred between the time they started runaway gas accretion and the time they 
contracted to their current dimensions and were able to gravitationally scatter planetesimals out of 
the Solar System. Thus, reasonable values of Mz for Jupiter and Saturn are ~ 10-30 Me and 
10-20 M© t respectively. In the cases of Uranus and Neptune, reasonable values for Mz would be 
somewhat less than their current high-Z masses at a time when the low-Z mass (hydrogen and 
helium) falls in the range 1-2 M© A reasonable value of Mz for these two planets is ~ 10 M© 

The results suggest that the solar nebula dissipated while Uranus and Neptune were in the second 
phase, during which, for a relatively long time, the masses of their gaseous envelopes were small 
but not negligible compared to the total masses. Our estimates for the formation time for Uranus 
fall in the range of 2-16 million years, while those for Jupiter and Saturn are 2—10 million years. 
The wide range corresponds to the results of an extensive series of tests in which key parameters 
were varied. Note that our estimates of formation times and surface densities follow from the 
properties of our Solar System and do not necessarily apply to giant planets in other planetary 
systems. 

Our paper on this topic has been published by Icarus. A copy of this paper (Pollack et al. 1996) is 
included with this report. 


Progress 

Preliminary work was performed in order to clearly define feasible goals and procedures for our 
continued improvement of our planetary evolution code. We began making modifications to the 
code that include updating tables and boundary conditions as well as improving the method for 
mass deposition of the incoming planetesimals in the envelope. We considered the applicability of 
using our code to simulate an extra-solar planet. A detailed description will be presented here. 

Our working code for simulating giant planet growth has been quite successful in accounting for a 
number of basic planetary' properties. Nevertheless, the simplifying assumptions that were used in 
past simulations could have major impacts on our results. One of the improvements we made is 
to better account for the mass deposited into the envelope of the forming giant planet by the 
incoming planetesimals. A planetesimal dissolution code is used to evaluate the planet's effective 
capture radius and the energy deposition profile of accreted material. In accord with the 
properties of comet Halley, we picture the planetesimal as consisting of small bits of rock and 
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organic matter embedded in a matrix of water ice (e.g., Jessberger et al. 1989). The ice acts as 
the “glue” which holds the planetesimal together. The surface temperature of the planetesimal is 
computed under the assumption of balance between heating and cooling, where heating includes 
gas drag and thermal radiation from the environment and cooling includes radiation emitted from 
the planetesimal surface and latent heat required to vaporize water ice. Vaporization occurs once 
the surface temperature exceeds a minimum vaporization temperature set by the vapor pressure 
(Podolak et al. 1988). When a layer of ice is vaporized, any rock or organics contained in that 
layer are also released into the envelope (referred to as “ablated material”). The fate of this 
ablated material then depends on the local ambient temperature, T env . When T env exceeds the 
vaporization temperature of the ice, rock, or organics, Tice > Track , and 7c//0.v> respectively), 
these materials vaporize, extracting energy from the layer in the case of rock and ice, and 
releasing energy in the case of organics. Otherwise, solid material keeps sinking slowly into the 
deeper regions of the envelope, releasing gravitational energy (through drag heating). Energy is 
also added to each mass shell corresponding to the conversion of kinetic energy into heat by the 
gas drag on the remaining planetesimal. Finally, the planetesimal is assumed to be fragmented 
into small (digestible) pieces when the gas dynamical pressure exceeds the compressional strength 
of the planetesimal. 

In our previous simulations the energy distribution for each mass shell of the envelope was 
calculated, but the dissolved mass was assumed to be deposited totally onto the core. The new 
modification enables us to compute models such that the dissolved mass will be self-consistently 
redistributed in the envelope. Our calculation will allow us to estimate the mass of the dissolved 
planetesimal material that remains in the envelope. The enrichment of the giant planet atmosphere 
in heavy elements with respect to the Sun will then be used to provide an additional check on 
models. 

The equation of state has been updated to that of Saumon et al. (1995). These tables were 
intended for application to low-mass stars, brown dwarfs, and giant planets and are for hydrogen 
and for helium in which nonideal effects are carefully included. In particular, pressure ionization 
of hydrogen has been explicitly treated. The opacity tables were updated to the new results of 
Alexander and Ferguson (1994). These Rosseland mean opacities have been computed over the 
temperature range 12,500—700 K for a range of different chemical compositions. The effects of 
atomic, molecular, and solid particulate absorbers and scatterers are included. 

Computation of the models described in Pollack et al. (1996) were stopped at the onset of 
runaway gas accretion. Improvements were made to the code in the boundary conditions to allow 
for hydrodynamic inflow of gas and to handle the late stages of evolution when the planet evolves 
at constant mass. A few runs were recently computed that simulate the termination of accretion 
of gas by the protoplanet (e.g. gap formation) and that follow the evolution through the final 
contraction and cooling phases, on time scales of 10^. The results have been compared with the 
models of Saumon et al. (1996) in the late stages of evolution, and there is good agreement. 
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New numerical simulations of the formation of the giant 
planets are presented, in which for the first time both the gas and 
pianetesimal accretion rates are calculated in a self-consistent, 
interactive fashion. The simulations combine three elements: 
(1) three-body accretion cross sections of solids onto an isolated 
planetary embryo, (2) a stellar evolution code for the planet’s 
gaseous envelope, and (3) a pianetesimal dissolution code 
within the envelope, used to evaluate the planet’s effective 
capture radius and the energy deposition profile of accreted 
material. Major assumptions include: The planet is embedded 
in a disk of gas and small planetesimals with locally uniform 
initial surface mass density, and planetesimals are not allowed 
to migrate into or out of the planet’s feeding zone. 

All simulations are characterized by three major phases. Dur- 
ing the first phase, the planet’s mass consists primarily of solid 
material. The pianetesimal accretion rate, which dominates 
that of gas, rapidly increases owing to runaway accretion, then 
decreases as the planet’s feeding zone is depleted. During the 
second phase, both solid and gas accretion rates are small 
and nearly independent of time. The third phase, marked by 
runaway gas accretion, starts when the solid and gas masses 
are about equal. It is engendered by a strong positive feedback 
on the gas accretion rates, driven by the rapid contraction of 
the gaseous envelope and the rapid expansion of the outer 
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boundary, which depends on the planet’s total mass. The overall 
evolutionary time scale is generally determined by the length 
of the second phase. 

The actual rates at which the giant planets accreted small 
planetesimals is probably intermediate between the constant 
rates assumed in most previous studies and the highly variable 
rates used here. Within the context of the adopted model of 
pianetesimal accretion, the joint constraints of the time scale 
for dissipation of the solar nebula and the current high-Z masses 
of the giant planets lead to estimates of the initial surface 
density ( aw) of planetesimals in the outer region of the solar 
nebula. The results show that a™ 10 g cm' 2 near Jupiter’s 
orbit and that a^ * cT l > where a is the distance from the Sun. 
These values are a factor of 3 to 4 times as high as that of 
the “minimum-mass” solar nebula at Jupiter’s distance and a 
factor of 2 to 3 times as high at Saturn’s distance. The estimates 
for the formation time of Jupiter and Saturn are 1 to 10 million 
years, whereas those for Uranus fall in the range 2 to 16 million 
years. These estimates follow from the properties of our Solar 
System and do not necessarily apply to giant planets in other 

planetary systems. © 1996 Academic Press. Inc. 


1. INTRODUCTION 

Unlike the terrestrial planets, the giant planets formed 
from significant quantities of both the gas and the solid 
material of the solar nebula. However, the giant planets 
are not made of solar proportions of the elements. Rather, 
all four giant planets preferentially accreted refractory ma- 


62 

0019-1035/96 $18.00 

Copyright © 1996 by Academic Pres s, Inc. 

All rights of reproduction in any form reserved. 




FORMATION OF GIANT PLANETS 


63 


terials, with the degree of enhancement, with respect to 
the Sun, varying progressively from a factor on the order 
of 5 for Jupiter to about 25 for Saturn to very roughly 300 
for Uranus and Neptune (e.g.. Pollack and Bodenheimer 
1989, Podolak et al 1993). Thus, it seems likely that the 
formation of the giant planets involved the “binary” accre- 
tion of solid planetesimals. the same process by which 
the terrestrial planets formed (Safronov 1969). However, 
unlike the terrestrial planets, the giant planets grew mas- 
sive enough to capture large quantities of gas from the 
solar nebula. 

Mizuno etai (1978) and Mixuno (1980) were the first to 
show that the above conceptual model was able to account 
approximately for the relative amounts of high- and low- 
Z materials in the giant planets. We refer to the gaseous 
component, which is primarily H : and He, as the “low-Z” 
material, where Z indicates the atomic number. We refer 
to the solid material, which includes “rock,” “CHON,” 
and "ice,” as the “high-Z” material, even though “ice” and 
“CHON” include significant amounts of H. In particular, 
Mizuno and collaborators constructed a series of equilib- 
rium model planets having solid (high-Z) cores and gas- 
eous envelopes that joined smoothly with the solar nebula 
at their tidal radii. The models had low-Z envelopes that 
grew exponentially with increasing core mass, the two 
masses becoming approximately equal when the core mass 
reached a “critical value” of M- nt ~ 10 Me. Mizuno was 
unable to construct equilibrium models for larger envelope 
masses, which led him to suggest that the giant planets 
underwent a hydrodynamicai collapse when their core 
masses exceeded M cnt . Gas would have accreted very rap- 
idly during this phase. The value of M cnt was found to 
depend very insensitively on the nebula boundary' condi- 
tions and weakly on the amount of grain opacity assumed 
to be present in the outer portions of the envelope. Thus, 
this “core instability” model appeared capable of ex- 
plaining why the high-Z masses of the giant planets were 
rather similar and had values on the order of 10 to 30Me- 
Bodenheimer and Pollack ( 1986, hereafter referred to 
as BP86) carried out the first evolutionary calculation of 
the core instability model. They constructed sequences of 
quasi-hydrostatic models that were connected in time by 
a prescribed rate of solid-body accretion and whose enve- 
lopes evolved in time due to gas accretion and radiation 
to space. They assumed that the solid-body accretion rate 
was time invariant for a given evolutionary sequence, al- 
though this rate was varied among the different sequences. 
In these simulations, the rate of gas accretion exceeded 
the rate of planetesimal accretion by an amount that grew 
exponentially with time once the core mass was sufficiently 
massive. Since this mass is not precisely determined by the 
simulations (the transition is fairly rapid, but not abrupt) 
and it corresponds approximately to the point where the 
high-Z and the low-Z masses are equal, it will be referred 


to as the “crossover mass” (M CTOSS ) in the remainder of 
this paper. During this “runaway” gas accretion phase, the 
envelope did not undergo a hydrodynamic collapse as long 
as the solar nebula could supply gas rapidly enough to 
compensate for an increasingly rapid contraction of the 
outer envelope and an increasingly rapid expansion of the 
planet's sphere of influence. BP86 obtained a value of 
M ctom (which they referred to as the critical core mass 
in analogy with Mizuno) — about 10 to 30Me — somewhat 
larger than the values of M^ obtained by Mizuno (1980). 
The values of A/ cr oss were found to be very insensitive to 
the boundary conditions with the solar nebula, in accord 
with Mizuno s (1980) values for M cm ; to be even more 
insensitive to the amount of grain opacity assumed in the 
outer envelope than Mizuno had found (due to the inclu- 
sion of water vapor opacity): and to have a mild sensitivity 
to the core accretion rate, with larger core accretion rates 
leading to a larger M CT0SS . Thus, one could speculate that 
the modest variation in high-Z masses among the giant 
planets reflected variations in their rates of planetesimal ac- 
cretion. 

Pollack et al (1986) and Podolak et al (1988) examined 
the ability of accreted planetesimals to pass through the 
envelopes of the giant planets and reach the core intact. 
Using the model envelopes of BP86, they found that a 
combination of gas drag, evaporation, and dynamical pres- 
sure made it increasingly difficult for planetesimals to ar- 
rive intact at the core boundary once the envelope mass 
exceeds a few percent of M ~ . Thus, a significant fraction 
of the planetesimals accreted by the giant planets should 
have been dissolved in their envelopes, enriching them in 
high-Z elements. Such a scenario is able to account in 
an approximate way for the observed enhancement (with 
respect to solar values) of some high-Z elements in the 
atmospheres of the current giant planets and for the pro- 
gressive enrichment of these elements from Jupiter to Sa- 
turn to Uranus/Neptune (Podolak et al. 1988. Simonelli 
etai 1989). 

Thus, it might appear that the core instability model is 
capable of explaining both the bulk mass properties and 
the atmospheric compositions of the giant planets. How- 
ever, this model faces several potentially very' important 
problems. First, Uranus and Neptune do not fit nicely into 
this picture. Their low-Z masses equal only about 10 to 
20% of their high-Z masses. Thus, they would not have 
been expected to have attained M cros s before their accretion 
was halted. It is not obvious why their high-Z masses should 
be similar to those of Jupiter and Saturn. Furthermore, 
the period during which accreting giant planets have low- 
and high-Z masses similar to those of Uranus and Neptune 
represents only a tiny fraction of their total accretion time 
in the calculations of BP86. While one could postulate that 
the solar nebula vanished at just the right time to account 
for one of these planets, the a priori probability that planets 
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with the properties of Uranus and Neptune would be found 
in the same system would be incredibly small in this sce- 
nario. 

The core instability scenario has been further analyzed 
by Wuchterl (1991), who used a radiation-hydrodynamics 
code rather than a quasi-static one. He finds that once the 
envelope mass has become comparable to the core mass, 
a dynamical instability develops that results in the ejection 
of much of the envelope, leaving a planet with a low-mass 
envelope and properties similar to those of Uranus and 
Neptune. The model does not account for the formation 
of Jupiter and Saturn, a problem that has not yet been 
resolved. Tajima and Nakagawa (preprint) have reexam- 
ined the evolution with the same assumptions as those of 
BP86 but with an independent code and have used a linear 
stability analysis to examine the properties of the envelope 
at all times. They find that the envelope is dynamically 
stable for all masses up to that of Jupiter, and that therefore 
the quasi-static approximation is justified. Nevertheless, 
Wuchterl (1995) continues to find dynamical instability 
unless the solar nebular density. p neb , is increased to 10~ 9 
g cm' 3 or higher (an order of magnitude higher than the 
standard value at Jupiter’s distance), in which case he finds 
stable accretion to high envelope masses. 

Another possible problem with the core instability hy- 
pothesis is that it does not account for the observed parti- 
tioning of high-Z material between a truly segregated inner 
core and the envelope. Recent interior models of the giant 
planets suggest that the cores of Jupiter and Saturn contain 
only a few M9, with the vast majority of the high-Z material 
residing in the envelopes (Zharkov and Gudkova 1991, 
Chabrier ex al 1992). Consequently, previous calculations 
that have implicitly assumed that planetesimals reach the 
core intact (e.g., Mizuno 1980. BP86) may not be directly 
relevant for estimating the mass of solids that needs to 
be accreted before runaway gas accretion takes place. In 
particular, these models probably overestimate the energy 
released by planetesimal accretion and thereby artificially 
delay the onset of rapid envelope contraction. 

Finally, there is a possible problem with the accretion 
time scale. In his classical calculations of planetesimal ac- 
cretion. Safronov (1969) obtained accretion time scales for 
Neptune that exceeded the age of the Solar System. A 
previous approach to the time scale problem was that of 
Stevenson (1984: see also Lissauer etal 1995), who consid- 
ered a core of the mass of Ganymede, from which icy 
material evaporated, forming a dense H 2 0-H 2 envelope 
which had a relatively small value of M CTOSS . However, other 
calculations suggest that this problem may be alleviated 
by some or all of the following factors: (1) rapid ‘"runaway” 
accretion of solids by the largest planetary embryos (Levin 
1978, Greenberg ex al 1978); (2) the possibility that the 
mass density of planetesimals in the giant planet region of 
the solar nebula exceeded somewhat the values given by 


the so-called “minimum mass” solar nebula (Lissauer 
1987); (3) more rapid accretion times found in multiple 
zone simulations of planetesimal accretion by Wetherill 
(see. e.g., Lissauer et al. 1995). These points also suggest 
the real possibility that the rate of planetesimal accretion 
may have deviated by wide margins from a time -invariant 
value, especially in the case when a single dominant mass 
is present (Lissauer 1987). 

In this paper, we improve the core instability model by 
constructing evolutionary models of the formation of the 
giant planets that allow for the interactions of planetesi- 
mals with the envelopes of the giant planets, and in which 
the rate of planetesimal accretion is calculated rather than 
prescribed. To do the latter, we must, of necessity, choose 
a particular model of planetesimal accretion. We selected 
the accretion model of Lissauer (1987) for the following 
reasons. First, it offers a promising means of solving the 
time scale problems alluded to above. Second, it represents 
a contrasting extreme from the prior assumption of a con- 
stant planetesimal accretion rate and so offers an opportu- 
nity to examine whether there is a qualitative difference 
in the results and. it is hoped, to bracket reality. The calcu- 
lations carry the evolution beyond the crossover mass into 
the phase of rapid gas accretion: they do not include the 
final phase where gas accretion terminates; thus, they do 
not attempt to explain the final masses of Jupiter and 
Saturn. Preliminary reports of the results presented in this 
paper were given by Podoiak ex al (1993) and Lissauer 
etal (1995). 

2. PROCEDURE 

To simulate the concurrent gas and solid accretion of 
the giant planets, we used an evolutionary model having 
three major components: a calculation of the three-body 
accretion rate of a single dominant-mass protoplanet sur- 
rounded by a large number of planetesimals: a calculation 
of the interaction of accreted planetesimals with the gas- 
eous envelope of the growing giant protoplanet; and a 
calculation of the gas accretion rate using a sequence of 
quasi-hydrostatic models having a core/envelope structure. 
These three components of the calculation were updated 
every time step in a self-consistent fashion in which rele- 
vant information from one component was used in the 
other components. We now describe these components, the 
key input quantities, and the limitations of our simulations. 

2.1. Planetesimal Accretion 

The early growth of the terrestrial planets from a swarm 
of planetesimals is thought to have involved “runaway” 
growth, in which there was a single dominant mass that 
grew rapidly from the accretion of nearby planetesimals 
(Greenberg et al 1978, 1984). The time scale associated 
with this phase of the formation of the planets was short 
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because the relative velocity of the planetesimals, u rel , was 
small compared with the escape velocity from the surface 
of the embryonic planet, u esc , and hence the planet’s gravi- 
tational cross section far exceeded its geometrical cross 
section. The runaway ended when the embryonic planet 
substantially depleted its “feeding zone.” The final phase 
of accretion of the terrestrial planets involved interacting 
embryos, with much higher values of i> re i and hence much 
longer time scales. Lissauer (1987) suggested that the giant 
planets may have had sufficient material in their nearby 
feeding zones so that they were able to reach during 
the runaway planetesimal accretion phase. Our specific 
assumptions are (1) there is a single accreting planet with 
i>ret < y esc and (2) the surface density of the planetesimal 
disk is a few times as large as that of the “minimum mass” 
solar nebula. As a result of these assumptions, there is a 
runaway phase during which A/cross is not quite reached. 

According to the classical theory of Safronov (1969), a 
solid protoplanet grows by accreting planetesimals whose 
orbits cross its orbit at a rate given by 


dM p 

dt 


uR]crCLF g. 


( 1 ) 


In Eq. (1), A/p is the mass of the giant protoplanet, R c is 
its effective or capture radius, cris the surface mass density 
of planetesimals, fi is the orbital frequency, and F g is the 
ratio of the gravitational cross section to the geometric 
cross section (“gravitational enhancement factor”). Ac- 
cording to Kepler’s laws of motion, 

n sc a" 3 ' 2 , (2) 

where a is the semimajor axis of the protoplanet. In many 
models, the time scales for planet growth increase steeply 
with increasing distance from the Sun because of the de- 
pendence of fl on a and because cr is expected to decrease 
with increasing a (Safronov 1969, Weidenschilling 1977). 
Time scale estimates for the Uranus/Neptune region ex- 
ceed 10 9 years (Safronov 1969). However, time scales for 
giant planet formation that are consistent with the lifetime 
of the solar nebula may be achieved if cr is somewhat 
enhanced above its value for a so-called “minimum’’-mass. 
solar nebula and if core growth takes place preferentially 
during the runaway planetesimal accretion phase when F g 
can be very large (up to 10 4 ) (Lissauer 1987, 1993). 

We considered a forming giant planet surrounded by 
planetesimals, all having equal masses and radii, m p and 
r p , situated in a disk of spatially constant (but time varying) 
mass density a. It was assumed that m p < A/ p . Accurate 
values of F g for this situation have been obtained by 
Greenzweig and Lissauer (1992). They did so by per- 
forming a large number of three-body (Sun, protoplanet, 
and planetesimal) orbital integrations where the planetesi- 


mals had a Rayleigh distribution of eccentricities and incli- 
nations. In our simulations, we used analytical expressions 
for F g that they derived as fits to their numerical calcula- 
tions (cf. Greenzweig 1991): F g is a function of in, and 
d c , where i H is the rms value of the planetesimals* orbital 
inclinations, e H is the rms value of the planetesimals* orbital 
eccentricities, and d z is the planetary embryo’s effective 
capture radius for planetesimals. All three parameters are 
expressed in Hill sphere units 


d c = 


a . 
— i, 

(3a) 

Rh 

a 

Rh 6 ' 

(3b) 


(3c) 

Rh 


where R H , the Hill sphere radius, is given by 



where A/ 0 is the mass of the Sun. 

According to the calculations of Greenzweig and Lis- 
sauer (1990, 1992), the planetesimals’ inclinations are con- 
trolled by their mutual gravitational scatterings, and their 
eccentricities are determined by a combination of these 
scatterings and gravitational interactions with the proto- 
planet at distances comparable to its Hill sphere radius. 
(Note that because of the coherence of the orbital motions 
of the planetesimals and protoplanet, gravitational interac- 
tions well within the protoplanet’s Hill sphere are not effec- 
tive in pumping up the planetesimals’ eccentricities to large 
values.) In particular, we use the prescription 


t’escp 

' H “ v 3^h 
e H = max(2i H ,2) 


( 5 ) 

( 6 ) 


where u exp is the escape velocity from the surface of a plan- 
etesimal. 

The planet’s accretion (feeding) zone was assumed to 
be an annulus that extended a radial distance, a fl on either 
side of its orbit. According to the simulations of 
Greenzweig and Lissauer (1990; cf. Kary and Lissauer 
1994), a ( is well approximated by 


Of = V 12 t Rh . (2) 

Thus, the accretion zone grows as the planet gains mass 
(independent of whether the accreted mass is solid or gas). 
The mass of planetesimals in the accretion zone is assumed 
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to equal the initial mass of planetesimals in the (current) 
accretion zone minus the amount that has already been 
accreted by the protoplanet; radial migration of planetesi- 
mals into and out of the accretion zone is therefore ne- 
glected. Random scatterings are assumed to spread the 
unaccreted planetesimals within the protoplanet’s reach 
uniformly over its accretion zone, so that the formulas for 
£ g given in Appendix B of Greenzweig and Lissauer (1992) 
are applicable. 

2.2. Interaction of Planetesimals with the Protoplanet 

The presence of a gaseous envelope around the high-Z 
core of a forming giant plane: can enhance the capture 
radius R c , and can lead to the deposition of mass and 
energy within the envelope when its mass is large enough. 
More precisely, these effects begin to occur when an incom- 
ing planetesimal intercepts a mass of gas comparable to its 
own mass (Pollack et ai 1986). We used the orbit trajectory 
code of Podolak et ai (1988) to evaluate these types of 
interactions of planetesimals with the protoplanet’s enve- 
lope. Here we summarize the protocols used. 

According to the calculations of Safronov (1969), the 
initial relative velocity of a planetesimal far from the pro- 
toplanet, on average, is given by 

u, = V T~* Vk ’ (8) 

where v k is the protoplanet's Kepierian velocity about the 
Sun. This velocity was divided by V2 to approximately 
account for the greater accretion rates of those planetesi- 
mals in the velocity distribution that have lower e and 
i. For various trial values of the impact parameter, the 
trajectory program used the analytical solution of the two- 
body problem (planetesimal, planet) for no gas drag to 
determine the planetesimal’s velocity vector at the point 
at which it reached the outer boundary of the protoplanet, 
R pt which is evaluated in the protoplanet structure code 
(cf. Subsection 2.3). Once the planetesimal entered the 
protoplanet's envelope, its trajectory was found by a nu- 
merical integration of the equations of motion with allow- 
ance for the gravitational field of the core and envelope 
and for gas drag. We used the dependence of the gas drag 
force on Mach and Reynolds numbers that is given in 
Podolak et ai (1988). The equations of motion were inte- 
grated with a fourth-order Runge-Kutta scheme, whose 
time step equalled a fraction of the local Kepierian period. 
This choice of time step ensured that smaller time steps 
were used close to the core, where more temporal resolu- 
tion was desirable. 

Critical values for the impact parameter and the associ- 
ated value of R c (= periapsis altitude) were obtained in 
an iterative fashion by finding the largest value of the 


impact parameter for which the planetesimal was captured. 
The criterion for capture was that at the end of the plane- 
tesimal’s first pass through the protoplanet’s envelope, its 
total energy (kinetic plus gravitational) was less than a 
small negative number. This number, £ esc , was set by the 
condition that the planetesimal had enough energy at R H 
to escape into a solar orbit: 

£ csc = -3m p CL 2 Rl i . (9) 

Strictly speaking, this equation applies to escape along 
the Sun-protoplanet line. However, the minimum energy 
needed for escape is close to zero in ail directions. 

Having established the value of the critical impact pa- 
rameter, we next carried out a set of trajectory calculations 
for a series of impact parameters that lay between zero 
and the critical value. For each choice of impact parameter, 
we followed the trajectory of the planetesimal until it either 
reached the surface of the core or totally vaporized in the 
protopianet’s envelope. Over the course of the trajectory, 
we kept track of the amount of mass vaporized within 
each mass shell of the envelope and the amount of energy 
deposited into each mass shell. We then averaged these 
results over the ensemble of impact parameters to deter- 
mine a mean value for mass and energy deposition in each 
shell of the protoplanet and at the core interface. 

In detail, our protocol for evaluating the mass and energy 
deposition profile of a planetesimal along its trajectory 
through the protopianet’s envelope was as follows. In ac- 
cord with the properties of comet Halley, we pictured the 
planetesimal as consisting of small bits of rock and organic 
matter embedded in a matrix of water ice (e.g., Jessberger 
et ai 1989). In this case, the ice acted as the “glue” that 
held the planetesimal together. The surface temperature 
of the planetesimal was computed under the assumption 
of balance between heating and cooling, where heating 
includes gas drag and thermal radiation from the environ- 
ment and cooling includes radiation emitted from the plan- 
etesimal surface and latent heat required to vaporize water 
ice. The emissivity of the small grains was assumed to be 
unity. Vaporization occurred at a rate set by the surface 
temperature and the associated vapor pressure (Podolak 
et at. 1988). When a layer of ice was vaporized, any rock 
or organics contained in that layer were also released into 
the envelope (referred to as “ablated material”). Their 
fate then depended on the local ambient temperature, T env . 
When Tenv exceeded the vaporization temperature of the 
ice, rock, or organics, ( T ice , r rock , and Tchon , respectively), 
these materials vaporized, extracting energy from the layer 
in the case of rock and ice, and releasing energy in the 
case of organics (as a result of chemical reactions with the 
ambient gas). Otherwise solid material kept sinking slowly 
into the deeper regions of the envelope, releasing gravita- 
tional energy (through drag heating). Energy was also 
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added to each mass shell corresponding to the conversion 
of kinetic energy into heat by the gas drag on the remaining 
planetesimal. Finally, the planetesimal was assumed to be 
fragmented into small (digestable) pieces when the gas 
dynamical pressure exceeded the compressional strength 
of the planetesimal. 

Based on the above discussion, the amount of energy 
released into a given mass shell i by the passage of a 
planetesimal and its associated debris, A £,, is given by 


A£j = £ d dsi + 2 Xj A/ni(0.5L*5 - L ; ^) 

'" l (10) 


3 i 




1 1 


where E d is the drag force exerted on the planetesimal in 
layer L ds t is the path length through layer i, Xj is the mass 
fraction of planetesimal constituent /, A m f is the total mass 
of the planetesimal vaporized and ablated in shell z\ u p is 
the local velocity of the planetesimal, L } is the latent heat 
of phase change of constituent ;, S :j is a Kronecker delta 
that equals 1 when constituent j undergoes a phase change 
in layer i and is 0 otherwise. A m.- is the total mass of the 
planetesimal ablated in shell f . <5- ; is the Kronecker delta 
that equals 1 when constituent ; is ablated in layer i T and 
reaches layer /, G is the gravitational constant, is the 
mass interior to mass shell i (envelope plus core), R t is the 
distance of mass shell i from the protoplanet’s center, and 
S"j is the Kronecker delta that equals 1 when constituent 
j ablates in mass shell i f and vaporizes in mass shell L The 
first term on the right-hand side of Eq. (10) represents 
heating of a mass shell by gas drag slowing of the planetesi- 
mal, the second term represents heating due to the dissipa- 
tion of the kinetic energy of the ablated material and cool- 
ing due to phase changes, and the third term represents 
heating due to gravitational energy release by sinking ma- 
terial that has ablated in layers above layer i and cooling 
by phase changes of this material. 

In applying the above equation for the energy added to 
each mass shell of the envelope, we considered two limiting 
cases concerning the ultimate fate of the ablated material. 
In fact the result will depend on whether the shell is radia- 
tive or convective and on the relative time scales of settling 
and mixing. On the one hand, once a given constituent of 
the ablated material is vaporized, it may be rapidly mixed 
with the surrounding H- and He-rich gas, in which case 
it will remain within the mass shell where it vaporized. 
Alternatively, mixing with environmental gas may be suf- 
ficiently sluggish that vaporized material continues to sink 
because its molecular weight is greater than that of the 
envelope. In this case, we allow vaporized material to sink 
all the way to the core interface, releasing gravitational 
energy as it does so. For brevity, we refer to these cases 


as the “no sinking” and “sinking cases. In either case, 
we add the mass of the planetesimal to the core, since our 
protoplanet structure code is not yet equipped to handle 
compositionally varying equations of state and opacity. 
Thus, in this one respect, the no sinking case is not strictly 
self-consistent. The effect of the dissolved material on the 
overall structure could be significant up to the time of 
crossover, and it will be considered in future calculations. 

Any remnant planetesimal that intersects the core re- 
leases its kinetic energy as heat at the core’s interface and 
uses up energy in phase changes that involve latent heat. 
We smear the net heating from this source and sink over 
a distance of one core radius into the envelope for reasons 
of numerical stability, as was done in our earlier calcula- 
tion (BP86). 

Table I summarizes the physical and chemical properties 
of the planetesimals used in our calculations. They are 
based on the most common types of materials found in 
comets, with special emphasis on the in situ measurements 
of comet Halley by the Giotto and Vega spacecraft (Jess- 
berger et ai 1989, Pollack et ai 1994). It is fortunate that 
the results of this paper do not depend sensitively on the 
precise properties given in this table, given that the compo- 
sition of the average comet is not known and may not 
represent a precise analog of the high-Z material from 
which the planets formed. 

2.3. Gas Accretion 

We constructed a time series of quasi-equilibrium core/ 
envelope models of forming giant planets to determine the 
rate at which gas was accreted from the surrounding solar 
nebula. The mass and radius of the core were set by the 
cumulative mass of planetesimals that had been accreted 
up until the time of current interest and by the assumed 
density of the core, p con . A value of 3.2 g/cm 3 was used 
for pcore, in accord with the materials composing the plane- 
tesimals and the high pressures and temperatures at the 
core interface (BP86). Our results do not depend sensi- 
tively on this choice. 

We used the same set of equations of state and opacity 
coefficients for the envelope gases as were used in BP86. 
The equations of state allow for dissociation, ionization 
(including H metalization), and nonideal gas effects and 
are based on detailed thermodynamical calculations (Gra- 
boske et ai 1975, Grossman et ai 1980). These equations 
of state apply to a solar mixture of elements, X = 0.74, 
Y = 0.243, Z = 0.017. Our opacity sources included small 
grains made of water ice, silicates, and iron for tempera- 
tures up to 1700 K, molecules (H 2 0, TiO) for temperatures 
up to 3000 K, and normal stellar sources at still higher 
temperatures (Alexander 1975, Alexander et ai 1983, Cox 
and Stewart 1970). These opacities are based on a solar 
mixture of elements; in particular, a solar abundance of 
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grains with approximately an interstellar size distribution 
is assumed. They do not include effects of organic grains 
or pressure-induced transitions of molecular hydrogen. 
Updated opacities for the range 800-10,000 K (Alexander 
and Ferguson 1994) were not yet available at the time 
these calculations were made, but these improved opacity 
estimates will be included in future calculations. 

Since we might expect that most of the small grains 
initially present in the outer part of the solar nebula would 
have already accreted into larger objects at the time of the 
start of our calculations, it may appear to be inconsistent 
for us to use the opacity of a solar mixture of small (less 
than a few tens of microns) grains. However, we point out 
that a significant fraction of the mass of the grains in the 
coma of comets have sizes in the "smair size regime (Ma- 
zetsera/. 1987, McDonnell etal. 1987). Thus, large amounts 
of small grains should have been released by planetesimal 
ablation in the outer envelopes. Furthermore, collisions of 
planetesimals produced ejecta containing small panicles 
in the surrounding solar nebula. Clearly, however, it is 
difficult to estimate the amount of small grains present 
in the outer envelopes of the forming giant planets (see 
Lissauer et ai 1995 for a more detailed discussion of this 
important topic). Fortunately, the simulations of giant 
planet formation by BP86 suggest that key results, such as 
the value of the M cross , do not depend sensitively on the 
amount of grain opacity. In particular, reducing the grain 
opacity by a factor of 50 led to only a 25% reduction in 
Mcross ; a further test reported below confirms the insensitiv- 
ity of Mcross but shows that the evolutionary time scale can 
be strongly affected. Finally, we note that grain opacity 
exceeds that due to H : as long as the abundance of small 
grains is more than 10~ 4 that for a solar mixture and that 
including the opacity of organic grains would boost the 
grain opacity at low and intermediate temperatures (<650 
K) by about a factor of 2 (Pollack et ai 1994). 

Quasi-equilibrium models of the envelope were con- 
structed by using the conventional stellar structure equa- 
tions of mass and energy conservation, hydrostatic equilib- 
rium, and the diffusion equation for radiative transfer 
(BP86). In convection zones the temperature gradient was 
approximated by the adiabatic gradient. The energy equa- 
tion includes three sources: the heat generated by captured 
planetesimals, PdV work from compression by gravita- 
tional forces, and cooling from release of internal heat. 
The boundary conditions at the inner edge of the envelope 
are that the luminosity is zero, the mass equals the core 
mass, and the radius equals the core radius. The outer 
radius of the envelope, R p , is set equal to the smaller of 
the Hill sphere (tidal) radius, R H . defined in Eq. (4), and 
the accretion radius, /? a , given by 


where c is the sound speed in the solar nebula. As discussed 
in BP86, a gas parcel located outside of has more ther- 
mal energy than gravitational energy binding it to the pro- 
toplanet. Hence, it is not part of the planet beyond /? a 
(or more precisely, it is no longer appropriate to use the 
equation of hydrostatic equilibrium). At R p , we require 
the envelope’s density and temperature to equal those in 
the surrounding solar nebula, p neb and T ncb . In actuality, 
this condition will not be met precisely at either R A or /? H » 
but rather will reflect the complicated flow of the solar 
nebula near a protoplanet. Fortunately, our results depend 
very insensitively on these outer boundary conditions (Mi- 
zuno 1980, BP86). 

Gas accretion occurs as a result of the contraction of 
the outer envelope and the steady increase in R p as the 
planet’s total mass increases. Gas from the surrounding 
solar nebula is assumed to flow freely into the evacuated 
volume at whatever rate is needed to restore the outer 
boundary conditions. Suppose that at time t the envelope 
has a radius R p (t) that is consistent with the outer boundary 
conditions. During time step A r. a planetesimal mass A m p 
is added to the planet, increasing the outer radius to R bd 
(where R bd = min[/? a , /? H ])« while the planet’s radius con- 
tracts to R p (t + A t). Thus, an amount of gas, Am ncb , will 
be added that is given by 

A/n„ eb = ^TrR 2 baPr^[ R bd - Rp{t + Af)]. (12) 

The mass of the added gas causes R bd to increase, and the 
gas is redistributed over the evacuated space in accord 
with the equation of hydrostatic equilibrium. Thus, we 
iteratively adjust the outer boundary and amount of added 
mass to obtain a self-consistent structure at time t + A t. 

2 A. Putting the Pieces Together 

We have now described the protocols used for the three 
main parts of the calculation. Here, we indicate how these 
pieces interact and summarize the key parameters used in 
our simulations. We begin the calculation at t = 0, with 
an initial model of the planet those total mass (almost 
entirely high-Z) is comparable to that of Mars. Using the 
procedures outlined in BP86. we find a quasi-equilibrium 
structure for the envelope of this initial model that matches 
the specified time-invariant values of p neb and r ncb at its 
outer boundary. We also specify the initial column mass 
density of planetesimals, cr^. their composition, and their 
radius, r p ; planetesimal radii and composition are assumed 
to remain constant over the course of accretion (an obvious 
oversimplication). At t = 0, the surface density of the disk 
is assumed to be constant and, thus, we do not take into 
account the small decrease in the planetesimal surface den- 
sity in the embryo feeding zone, which has already occurred 
as a result of the mass incorporated in the embryo. Tables 
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TABLE I 


Properties of Planetesimals 


Component® 

Property 

H 2 0 Ice 

Rock 

CHON 

Total 

mass fraction 

0.397 

0.308 

0.295 

1 

density (g/cm 3 ) 

0.92 

3.45 

1.5 

1.39 

latent heat 3 (erg/g) 

c 2.8 x lO 10 

=8.08 x 10 l ° 

d - 7.0 x 10 10 

1.54 x 10 10 

vaporization temperature (K) 

165 

1500 

650 



- The three major components of the planetesimals are water ice, ferromagnesium silicates (“rock”), 

and organics i**CHON M ). ^ TT _ vT • . 

b The latent heats of ice and rock are endothermic, whereas that of the CHON is exothermic. 


c Podolak er aL ( 1988). 
d Estimated. 

I and II summarize the key parameters that are used in 
the calculation and their nominal set of values. 

We next wish to find a new equilibrium model at the 
end of time step A r. determining the rates of planetesimal 
and gas accretion in the process. To do this, we first 
use the formulas for three-body accretion to determine 
the mass of accreted planetesimals. In so doing, we use 
the properties of the model of the planet at the beginning 
of the time step to determine R c and e H . Then, we use 
the trajectory code to evaluate the mass and energy 
deposited by the accreting planetesimals in each shell 
of the protoplanet's envelope and at the core interface. 
In doing this calculation, we use a preliminary' updated 
envelope structure that takes account of the added plane- 
tesimal and gas masses (and an associated rezoning of 
the mass shells). For a particular sequence, we choose 
either the sinking or no-sinking option for the vaporized 
material from the planetesimals in the envelope. Finally, 
the energy profile found from the trajectory calculations 
is used to calculate a final equilibrium structure for 
the envelope. 


TABLE n 

Key Model Parameters and Their Nominal Values 


Parameter 

Nominal Value 

orbital distance 

5.2 A. U. 

planetesimal radius 

100 km 

other planetesimal properties 

see Table I 

initial planetesimal surface density 

10 g/cm 2 

fate of dissolved planetesimal 

sinks to core interface 

nebula temperature 

150 K 

nebula density 

5.0 x 10~ u g/cm 3 


We now have a new model for the protoplanet and are 
ready to take the next time step. First, we readjust the 
column mass density of planetesimals to allow for the mass 
that has been accreted in the previous time step and the 


TABLE III 
Input Parameters 


case 

(g/«m 2 ) 

&\mt,XY 

(g/cm J ) 

r ? 

(km) 

a 

(A.U.) 

Tnzo 

(K) 

Pnei 

(g/«n 3 ) 

i. 

J1 

10. 

700. 

LOO 

5.203 

150 

5.0 x L0 -11 

1 

Jla* 

10. 

700. 

100 

5.203 

150 

5.0 x 10 -u 

1 

Jib 6 

10. 

700. 

100 

5.203 

150 

5.0 x 10 _u 

1 

JlC c 

10. 

700. 

100 

5.203 

150 

5.0 x L0“ n 

1 

J2 

7.5 

525. 

100 

5.203 

150 

5.0 x 10 _u 

1 

J3 

15. 

1050. 

100 

5.203 

150 

5.0 x 10 _u 

1 

J4 

10. 

700. 

100 

5.203 

150 

5.0 x 10 -u 

0 

J5* 

10. 

700. 

100 

5.203 

150 

5.0 x 10 -u 

1 

J6* 

10. 

700. 

100 

5.203 

150 

5.0 x 10 -u 

1 

J7 

10. 

700. 

1 

5.203 

150 

5.0 x 10“ n 

1 

J8 

10. 

700. 

1 

5.203 

150 

5.0 x 10 -11 

0 

SI 

3. 

210. 

100 

9.539 

100 

2.5 x 10 _u 

1 

S2 

3. 

210. 

100 

9.539 

100 

2.5 x 10 -u 

0 

U1 

0.75 

52.5 

100 

19.18 

75 

1.0 x 10 _u 

1 

TJ2 

0.75 

52.5 

1 

19.18 

75 

1.0 x 10" u 

1 


•Planetesimal accretion arbitrarily stopped at L5 myr. 
b Planetesimal accretion arbitrarily stopped at 3.5 myr. 
c Planetesimal accretion arbitrarily stopped at 6.8 myr. 

4 All planetesimals reach the core and energy is deposited within 
one core radius of core-envelope boundary. 

• Grain opacity in the envelope equals 2% of nominal (solar compo- 
sition) value. 
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TABLE IV 
Results 


FIRST MAXIMUM END OF PHASE 1 CROSS-OVER POINT ENDPOINT 


ran 

Time* 

Mxy" 

M z 

Tune 

Mxy 

Mz 

Mxy=Mz c 

Time 

Mcrota 

M X y 

M z 

Time 

Mxy 

M z 

Mxy 

Mz 

J1 

0.48S 

3.3xl0“ 3 

6.96 

0.01 

0.16 

11.4 

5.1x10-* 

7.58 

16.17 

1.2xl0“ 5 

2.7x10“* 

8.00 

64.4 

21.5 

7.3xl0“ 3 

3.7xl0“ 4 

Jl« 

tt 

tt 

it 

it 

tt 

tt 

// 

3.32 

12.24 

3.0x10“* 

0 . 

3.54 

29.8 

12.2 

3.3xI0“ 4 

0 . 

Jib 

it 

tt 

it 

it 

tt 

n 

// 

4.60 

13.04 

6.1xl0“ 5 

0 . 

4.75 

34.3 

13.0 

l.OxlO” 3 

0. 

Jlc 

it 

tt 

it 

it 

tt 

it 

It 

7.07 

14.94 

9.7xl0“ 5 

0 . 

7.16 

34.3 

14.9 

2.7xl0" 3 

0. 

J2 

0.580 

8.4x10“* 

4.65 

0.81 

0.08 

7.4 

1.0x10“* 

48.0 

10.51 

1.3x10"* 

2.3xlO“ T 

50.0 

15.8 

11.4 

1.1x10“* 

6.9x10“* 

J3 

0.388 

1.6xl0“ 1 2 

12.1 

0.45 

0.65 

21.0 

2.5x10"* 

1.51 

29.61 

L2xl0“ 4 

3.1x10“* 

1.57 

56.0 

33.8 

1.1x10“ 3 

1.5xl0- 4 

J4 

0.485 

9.1xl0" 3 

7.2 

0.59 

0.30 

11.5 

6.9x10“* 

1.86 

10.07 

0.3x10“* 

1.8x10“* 

1.93 

48.2 

19.7 

5. 1x10“ 3 

3.7xl0“ 4 

J5 

0.489 

4.0xl0- 3 

7.18 

0.63 

0.23 

11.5 

3.4x10"* 

7.65 

10.17 

1.2x10“* 

2.9x10“* 

8.00 

38.2 

19.3 

1.5xl0“ 3 

9.0x10“* 

J6 

0.405 

l.lxlO" 2 

7.56 

j 0.55 

0.30 

11.5 

8.7xl0“ 6 

2.75 

18.18 

1.8x10“* 

3.9x10“* 

2.88 

19.4 

16.8 

2.6x10“* 

5.2x10“* 

J7 

0.214 

7.9xl0” 4 

5.64 

i 0.26 

0.14 

11.5 

7.1x10“* 

8.94 

18.18 

5.0x10“* 

2.8x10“* 

7.31 

44.7 

20.22 

l.lxlO" 3 

l.lxlO" 4 

J8 

0.214 

3.0xlQ“ 3 

5.44 

j 0.28 

0.62 

11.6 

1.4xl0“ 5 

1.22 

10.07 

2.6xl0“ 4 

3.3x10"* 

1.25 

34.2 

18.6 

2.5xl0“ 3 

2.4xl0“ 4 

SI 

1.90 

7.9xl0' 3 

7.49 

: 2.48 

0.45 

11.7 

1.5x10“* 

9.50 

16.34 

LOxlO" 5 

2.7x10“* 

9.80 

23.2 

17.5 

4.2x10“* 

7.0x10“® 

S2 

1.87 

4.6xl0" 2 

7.31 

2.23 

1.11 

11.7 

5.0x10“* 

3.24 

15.87 

l.lxlO" 4 

1.7x10“* 

3.29 

27.0 

17.33 

6.8xl0“ 4 

5.3x10“* 

U1 

12.9 

3.4xl0" 2 

7.32 

15.2 

0.83 

11.7 

8.7x10“ t 

21.9 

16.24 

2.8x10“* 

2.6x10"* 

22.05 

21.9 

16.86 

8.3x10“* 

6.3x10“® 

U2 

0.476 

2.3xl0“ 3 

5.48 

0.92 

0.48 

11.9 

1.8x10“* 

6.69 

16.57 

1.2x10“* 

3.3x10“® 

6.94 

26.9 

18.07 

1.7xl0“ 4 

1.6x10“* 


fl Time is in units of millions of years, myr. 

b Mass is in units of Earth’s mass. .\/ e . 

c Accretion rate is in units of Earth masses per year, A/ 9 /year. 


expansion of the feeding zone (see Section 2.1). Then, we 
follow the same sequence of steps to obtain a new model 
of the protoplanet and the new rates of planetesimal and 
gas accretion. 

2.5. Key Assumptions 

Here, we summarize key assumptions made in our simu- 
lations and define their basic limitations. These include: 

1. The opacity in the outer envelope is determined by 
a solar mixture of small grains. We will comment below 
on the effects of a change in the abundance of small grains. 
We also assume solar abundances in calculating the opacity 
in deeper regions of the envelope where molecular opacit- 
ies dominate. Here, we may have underestimated the true 
opacity throughout much of the accretion (e.g., enhanced 
amounts of dissolved H 2 0 would raise the opacity). 

2. The equation of state for the envelope is that for 

a solar mixture of elements. This will start to become a 

questionable assumption when a large amount of planetesi- 
mal mass has dissolved in the envelope. For example, the 

composition gradient introduced by the distribution of dis- 
solved heavy material could affect the extent of convection 
zones. But, variations in the equation of state resulting 

from the addition of dissolved high-Z material probably 
affected the evolution of the planet far less than the 


changes in the opacity due to the same addition of high- 
Z material. 

3. During the entire period of growth of a giant planet, 
it is assumed to be the sole dominant mass in the region 
of its feeding zone, i.e., there are no competing embryos, 
and planetesimal sizes and random velocities remain small. 
A corollary of this assumption is that accretion can be 
described as a quasi-continuous process, as opposed to a 
discontinuous one involving the occasional accretion of a 
massive planetesimal. 

4. Planetesimals are assumed to be well mixed within 
the planet’s feeding zone, which grows as the planet’s mass 
increases, but planetesimals are not allowed to migrate 
into or out of the planet’s feeding zone as a consequence 
of their own motion. Tidal interaction (Lin and Papaloizou 
1993) between the protoplanet and the disk, or migration 
of the protoplanet, is not considered. It is not at all obvious 
that these various assumptions are valid, but no well- 
defined quantitatively justifiable alternative assumptions 
are available. 

5. Hydrodynamic effects are not considered in the evo- 
lution of the envelope. Although Wuchterl (1991) found 
that dynamical instability occurred in his models once the 
envelope mass became comparable to the core mass, most 
of the present results are based on the phases before the 
crossover mass is reached. 
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3. RESULTS 

In this section, we first present the results for a baseline 
simulation in some detail and then examine the sensitivity 
of the results to variations in key parameters. Table III 
states the parameters for each case and Table IV gives 
basic results. In all cases, the simulations begin with a 
protoplanet having a mass comparable to that of Mars 
with almost all its mass in a high-Z core. We continue the 
evolution past the onset of runaway gas accretion, which 
operationally is defined by the point where the gas accre- 
tion rate exceeds the planetesimal accretion rate by a factor 
of 10 and is increasing in a quasi-exponential fashion 
with time. 

We judge the applicability of a given simulation to plan- 
ets in our Solar System using two basic yardsticks. One 
yardstick is provided by the time required to reach the 
runaway gas accretion phase. This time interval should be 
less than the lifetime of the gas component of the solar 
nebula. r sn , for successful models of Jupiter and Saturn 
and greater than f sn for successful models of Uranus and 
Neptune. Limited observations of accretion disks around 
young stars suggest that r sn ^ 10 - years, based on observa- 
tions of the dust component. The lifetime of the gas compo- 
nent is less well constrained observationally (Strom et al 
1993). 

A second yardstick is provided by the amount of high- 
Z mass accreted, A/ z . In the case of Jupiter and Saturn, 
M z at the end of a successful simulation should be compa- 
rable to, but somewhat smaller than, the current high- 
Z masses of these planets, since additional accretion of 
planetesimals occurred between the time they started run- 
away gas accretion and the time they contracted to their 
current dimensions and were able to gravitationally scatter 
planetesimals out of the Solar System. Thus, reasonable 
values of M z for Jupiter and Saturn are 10-30 A/® and 
10-20 A/®, respectively. In the cases of Uranus and Nep- 
tune, reasonable values for M z would be somewhat less 
than their current high-Z masses at a time when the low- 
Z mass Mxy falls in the range 1-2AJ®. A reasonable value 
of M z for these two planets is —10 M©. 

3.1. Baseline Model for Jupiter 

Our baseline model (case Jl). which uses the parameters 
given in Tables I and II, is meant to provide a reasonable 
simulation of the formation of Jupiter, as judged by the 
yardsticks defined above. Figure la shows the evolutionary 
behavior of the masses Mxr and M z . (Note that values 
quoted for Mxy refer to ail of the mass accreted in the gas 
phase and, thus, they include a small fraction of heavy 
elements.) Figure lb shows the planetesimal accretion rate 
( dM z ldt ) and the gas accretion rate (dMxvIdt) for this 
baseline model, and Fig. 1c illustrates the luminosity. Ac- 
cording to these figures, there are three main phases to 


the accretion of our model Jupiter. Phase 1 is characterized 
by rapidly varying rates of planetesimal and gas accretion. 
Throughout phase 1, dM z ldt exceeds dMxyldt. Initially, 
there is a very large difference (many orders of magnitude) 
between these two rates; however, they become progres- 
sively more comparable as time advances. Over much of 
phase 1, dM z ldt increases steeply. After a maximum at 
5 x 10 s years, it declines sharply. Meanwhile, dMxyldt 
keeps steadily growing by many orders of magnitude from 
its extremely low initial value. 

The second phase of accretion is characterized by rela- 
tively time-invariant values of dM z l dt and dMxyldt , with 
dMxyldt > dM z ldt. We note that the small fluctuations in 
the accretion rates that are particularly noticeable during 
this phase (Fig. lb) are a numerical artifact that stems from 
the iterative scheme used to adjust R p (the actual planetary 
radius) to have a value approximately equal to R M . Finally, 
phase 3 is defined by rapidly increasing rates of gas and 
planetesimal accretion, with dMxyldt exceeding dM z ldt 
by steadily increasing amounts. 

Insights into the physics that controls the accretion rates 
during the three phases (but especially phase 1) may be 
obtained if one examines the evolutionary behavior of 
the surface density of planetesimals (Fig. Id) within the 
protoplanet’s feeding zone. cr z . Initially, before cr z de- 
creases significantly, dM z >'dt is expected to increase rapidly 
due to an increase in the capture radius of the growing 
protoplanet (e.g., Lissauer 1987, Wetherill and Stewart 
1989). The planetesimal capture radius, R c , initially simply 
equals R cor e and 

(13) 

at 

when gravitational focusing is taken into account in the 
two-body approximation. But later, when the envelope 
becomes sufficiently massive, planetesimals are captured 
by gas drag and, hence, R c exceeds R CO rc by progressively 
larger amounts (cf. Fig. le). As a result, dM z ldt increases 
even more rapidly with time. 

A decline in dM z ldt begins when the cumulative amount 
of accreted high-Z material becomes a significant fraction 
of the mass initially contained within the current bound- 
aries of the feeding zone. Such a depletion is inevitable 
within the context of our assumptions about the source of 
planetesimals, since the mass of material within the feeding 
zone is roughly proportional to R H [Eq- 00] and x 
Mp 3 [Eq. (4)]. The sharp decline of cr z with time in Fig. 
Id is the result of a progressive and ultimately nearly com- 
plete depletion of planetesimals in the protoplanet’s feed- 
ing zone brought about by prior accretion. Thus, phase 1 
of our evolutionary simulations denotes the period during 
which runaway planetesimal accretion occurs, and it ends 
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FIG. 1. (a) Planet’s mass as a function of time for our baseline model, case Jl. In this case, the planet is located at 5.2 AU. the initial surface 

density of the protoplanetary disk is 10 g. cm 2 , and pianetesimals that dissolve during their journey through the planet’s envelope are allowed to 
sink to the planet’s core; other parameters are listed in Table III. The solid line represents accumulated solid mass, the dotted line accumulated 
gas mass, and the dot-dashed line the planet’s total mass. The planet’s growth occurs in three fairly well-defined stages: During the first X 10 5 
years, the planet accumulates solids by rapid runaway accretion; this lt phase 1” ends when the planet has severely depleted its feeding zone of 
pianetesimals. The accretion rates of gas and solids are nearly constant with Mxy ** 2-3 Mz during most of the -~7 X 10 6 years’ duration of phase 

2. The planet’s growth accelerates toward the end of phase 2, and runaway accumulation of gas (and, to a lesser extent, solids) characterizes phase 

3. The simulation is stopped when accretion becomes so rapid that our model breaks down. The endpoint is thus an artifact of our technique and 
should not be interpreted as an estimate of the planet’s final mass, (b) Logarithm of the mass accretion rates of pianetesimals (solid line) and gas 
(dotted line) for case Jl. Note that the initial accretion rate of gas is extremely slow, but that its value increases rapidly during phase 1 and early 
phase 2. The small-scale structure which is particularly prominent during phase 2 is an artifact produced by our method of computation of the 
added gas mass from the solar nebula, (c) Luminosity of the protoplanet as a function of time for case JL Note the strong correlation between 
luminosity and accretion rate (cf. b). (d) Surface density of pianetesimals in the feeding zone as a function of time for case Jl. Pianetesimals become 
substantially depleted within the planet's accretion zone during the latter part of phase 1, and the local surface density of pianetesimals remains 
small throughout phase 2. (e) Four measures of the radius of the growing planetary embryo in case Jl. The solid curve shows the radius of the 
planet’s core, R c ore , assuming all accreted pianetesimals settle down to this core. The dashed curve represents the effective capture radius for 
pianetesimals 100 km in radius, R c . The dotted line shows the outer boundary of the gaseous envelope at the ‘"end” of a timestep, R p . The long- 
and short-dashed curve represents the planet’s accretion radius, R t . 


when the protoplanet has virtually emptied its feeding zone 
of pianetesimals. 

If this simulation had been done in a gas-free environ- 
ment, as might be appropriate for the formation of the 
terrestrial planets, then the next phase would have to in- 


volve interacting embryos for accretion to reach the desired 
culmination point (Lissauer 1987, Lissauer and Stewart 
1993). However, it is possible to carry our simulations of 
the formation of the giant planets to a reasonable endpoint 
without involving interacting embryos, because of the im- 
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pact of gas accretion on the subsequent evolution. In partic- 
ular, the gaseous envelope is massive enough at the end 
of phase 1 for its contraction to lead to further gas accre- 
tion, which augments the planet’s total mass. This results 
in a progressive increase of its Hill sphere radius, and 
hence, the size of its feeding zone keeps increasing, bring- 
ing new planetesimais within its sphere of influence. Thus, 
phase 2 involves a controlled interaction between the ac- 
cretion rates of gas and planetesimais. The nature of phase 
2 is of course dependent on assumption 3 of Section 2.5, 
namely, that there are no competing embryos in the neigh- 
borhood so that the supply of planetesimais is continuous. 
If this condition did not exist, then gas accretion during 
phase 2 would lead to merging of nearby embryos rather 
than smooth accretion of planetesimais. We will defer until 
Section 4.1 a discussion of the mechanisms that are respon- 
sible for the differences in accretion style during phases 2 
and 3. Here, we simply note that phase 3 is analogous 
to the classical “runaway" gas accretion phase of earlier 
calculations (e.g., BP86). However, phase 2 is an entirely 
new phase that was not present in previous simulations 
of the formation of the giant planets. It represents the 
transition phase between runaway accretion of solids and 
runaway accretion of gas. 

As can be seen from Figs, la and b, phase 1 lasts only 
about 6 X 10 5 years. This time scale, fp h i, can readily be 
estimated from the set of equations (Lissauer 1987) 


M lso = C^aro-nit) 

CzMlso 

fphl (m^PJL)’ 


3/2 


(14) 

(15) 


density of planetesimais in the feeding zone: and other 
variables have been defined earlier. We obtained a large 
value of M iso by selecting a cr, mI that was somewhat larger 
than that given by a so-called minimum mass solar nebula 
and a small value for r piU by considering a situation where 
F, is very large even at r = 0 due to the small random 
velocities of the planetesimais (Lissauer 1987). 

At the end of phase 1. M z 12 in agreement with 
Eq. (14), and M x y < 1A4- Subsequently, in phase 2, M z 
increases by another 4 and Mxy increases to a value 
essentially equal to that of M z ■ During phase 2. the surface 
density of the solids remains very small. The planetesimal 
accretion rate is thus essentially equal to the rate at which 
planetesimais enter the planet s accretion zone. The initial 
mass of planetesimais within the planet’s accretion zone 
is proportional to the one-third power of the planet s mass, 
so during phase 2 


M z + -V/y y x M\. (lb) 

Differentiating Eq. (16) with respect to time, we obtain 
the following relationship between the accretion rates of 
solids and gas during phase 2: 

M z ~ (?. -3 -j^J Mxy ■ (I 7 ) 

The numerical results (Fig. lb) are consistent with this 
expression. It follows from Eq. (16) that the crossover 
mass, Af eross , at which M z = Mxy , is given by 

Across ^ V2 M lso . (18) 


where M lso is the planet’s “isolation mass” (its mass after 
its feeding zone has been depleted), Cj. = 1.56 X 10 25 g if 
cr init is in units of g/cm 2 and a is in units of AU, and C 2 = 
8.126 with all quantities in cgs units; <x init is the initial surface 


The inequality is present in Eq. (18) because some plane- 
tesimais reside within the planet’s accretion zone when the 
crossover mass is reached. Note that the crossover masses 
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FIG. 2. (a) Cumulative mass of iow-Z material (dots), high-Z material (solid), and total (dot-dash) as a function of time for case J2. The 

duration of phase 2 is a factor of 7 longer than in case J1 as a result of a 259b drop in initial surface mass density of planetesimals. and is inconsistent 
with the best available estimate of the lifetime of gas within the solar nebula, t m s 10 7 years, (b) Cumulative masses (as in a) as a function of time 
for case J3. The increase in the surface density of planetesimals by 50% relative to case J1 leads to a much more rapid progression through phase 
2; however, the crossover mass. 30 Af$, may be larger than the amount of condensible material in Jupiter. 


found in our simulations are always within a few percent 
of V2 times the mass at the end of phase 1 (cf. Table 4), 
except in those runs where we modified our procedure by 
terminating planetesimal accretion during phase 2. 

Phase 2 lasts about 7 X 10 6 years. Thus, for our baseline 
model, the time required for the protoplanet to reach run- 
away gas accretion (phase 3) is determined almost solely 
by the duration of phase 2. This time scale is comparable 
to the estimated lifetime of the solar nebula, f sn . In addition, 
the mass of high-Z material accumulated by the end of 
our simulation, 21 .5M$, is comparable to or somewhat 
less than current estimates of \t z for present-day Jupiter 
(Zharkov and Gudkova 1991. Chabrier et al 1992). We 
conclude that our baseline model is consistent with the 
two basic yardsticks for judging the reasonableness of a 
simulation. Additional comparisons with observational 
and theoretical constraints will be made in Section 4.3. 

3.2. Other Models for Jupiter 

We now focus our attention on the effects of parameters 
that significantly influence the results of our simulations 
of the formation of the giant planets. These parameters 
include cr mK , <5 S , and a . The parameter S s equals 1 when 
we allow dissolved planetesimal-derived material to con- 
tinue to sink toward the core and release gravitational 
energy; it equals 0 when this material is not allowed to 
sink. Effects of opacity, planetesimal size, and the outer 
boundary condition are also considered. 

Figures 2a and b illustrate the great sensitivity of the 
accretion rates to variations in cr^ from its baseline value 
of 10 g/ cm 2 . Decreasing cr^. to 7.5 g/cm 2 (case J2) or 


increasing it to 15 g/cm 2 (case J3) greatly alters both the 
time it takes the protoplanet to reach the runaway gas 
accretion phase and the mass of high-Z material that it 
contains at this point. In particular, the interval of time 
from the start to the finish of our simulations. r end , equals 
5.0 X 10 7 years for the low-o^ case (J2), 8.0 X 10 6 years 
for the nominal-surface-density case (Jl), and 1.6 x 10 6 
years for the high-cr mil case (J3). Thus, varying cr mu by only 
a factor of 2 results in a factor of 30 variation in f <nd . Since 
r phl is relatively insensitive to cr imt , this effect is almost 
entirely determined by the duration of phase 2. The value 
of M z equals 11.4, 21.5, and 3 3.8 A/ 9 at t = r end for the low-, 
nominal-, and high-surface-density cases, respectively. 
Thus, varying cr inu by a factor of 2 produces a factor of 3 
spread in M z , as expected from Eq. (14). 

The great sensitivity of r 5nd and M z to a mil makes it 
possible to place very tight constraints on the actual value 
of cr ma of the solar nebula, within the context of our basic 
assumptions. In particular, a ^ has to lie within a few tens 
of percent of 10 g/cm 2 at Jupiter’s distance from the Sun 
for our simulations to be consistent with the two basic 
yardsticks for reasonableness of our results. If cr iml is s7.5 
g/cm 2 , then J end exceeds t SD . If <r init > 15 g/cm 2 , then the 
value of M z at f end exceeds the current high-Z mass of Ju- 
piter. 

We next examine the sensitivity of our results to the 
value of S s (case J4). Figure 3 shows the evolutionary his- 
tory of the masses when S s = 0 (no sinking). Comparing 
these results with those shown in Fig. la for S s = 1, we 
see that t end is shortened by about a factor of 4 when the 
vaporized material is not allowed to sink. However, M z at 
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FIG. 3. Cumulative masses las in Figs, la and 2) as a function of 
time for case J4. The duration of phase 2 is less than in case J1 because 
dissolved planetesimals do not sink, so less energy is available to support 
the planet’s envelope against gravitational contraction. 

t = fend for S s = 0 (19.7 A/.,) is essentially the same as its 
value for S s = 1 (21.5A/*). considering the fact that the 
latter case was evolved somewhat further. Thus, the bounds 
on acceptable values of cr imi for the no-sinking case remain 
the same on the high end and are slightly lower on the 
low end. in comparison with the corresponding bounds for 
the sinking case. In connection with the parameter S s , 
another test, case J5, was performed with the same parame- 
ters as in case J 1 but under the assumption that all planetes- 
imals reached the core and deposited their energy within 
one core radius of the core boundary; this procedure is 
that followed by BP86. The total energy released by plane- 
tesimal accretion is thus the same as in case Jl, but the 
distribution in radius is somewhat different. The evolution 
is essentially the same as that for case Jl. 

In case J6, the grain opaciry in the envelope is reduced 
by a factor of 50, although the molecular opacity, which 
dominates at temperatures above 1700 K, remains the 
same. The results for f phl and .Vf cross are hardly changed, 
but the time to transit phase 2 is only 2.2 X 10 6 years for 
case J6 as compared with 7.0 x 10 6 years for case Jl. Thus, 
the overall evolutionary time scale is reduced by almost a 
factor of 3, indicating that the details of the opacity are in 
fact significant. 

The effect of the outer boundary condition p neb is tested 
by reducing it by a factor 10. leaving all other parameters 
the same as in case Jl. As expected from previous numeri- 
cal (BP86) and analytical (Stevenson 1982, Lissauer et al 
1995) results there is very litile effect on the evolution. At 
Mxy ~ the value of is practically identical to 

that in case Jl, and the evolutionary time is a mere 1.7 X 
10 5 years longer. 


The effect of changing the assumed planetesimal size is 
considered in cases J7 and J8, which are calculated with 
r p = 1 km and with S s = 1 and 0, respectively. In both 
cases the isolation mass at the end of phase 1 is the same 
as in case Jl, but the time is reduced by about a factor of 
2, as a result of an enhanced gravitational focusing factor 
at early times. The time spent in phase 2 by case J7 is 6.7 X 
10 b years, practically the same as in case Jl. while case J8 
spent 1 x 10 6 years in this phase, slightly shorter than the 
corresponding time for the analogous case J4. Because 
most of the time is spent in phase 2, the effect of the 
planetesimal size on the evolution of Jupiter is small, how- 
ever, it is much more important in the case of Uranus, 
which is discussed in the following subsection. 

3 J. Models for Saturn and Uranus 

We now assess the impact of distance from the Sun on 
our results by considering values of a that correspond to 
the current orbital semimaior axes of the giant planets. We 
first increase a from 5.2 to 9.5 AU, the value appropriate for 
Saturn (case SI ). Since the current high-Z masses of Jupiter 
and Saturn are similar (Zharkov and Gudkova 1991. Cha- 
brier et al 1992), we want to pick a value for the isolation 
mass of Saturn that is comparable to the isolation mass of 
Jupiter for our nominal case. According to Eq. (14) for 
Mi S o, we therefore need to scale oi nit approximately as a . 

Figure 4 shows the evolutionary history of the masses 
for our nominal Saturn model (cr imt = 3 g/ cm 2 ). Not surpris- 
ingly, phase 1 for our nominal Saturn model lasts about 
four times longer than for our nominal Jupiter model; 
a* imt n is smaller by a factor of 8 but F $ is larger by a factor 
of 2 (cf. [Eq. (1)]. However. r end is only slightly larger for 



FIG. 4. Cumulative masses (as in Figs, la and 2) as a function of 
time for Saturn model case SI. The duration of phase 2 and the crossover 
mass iVfcoss are similar to case Jl. because the planetary isolation mass 
and the input physics are similar. 
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FIG. 5. (a) Cumulative masses (as in Figs, la and 2) as a function of time for Uranus case Ul. Note the relatively long period during which 

the model's mass and composition are similar to those of present-day Uranus, (b) Cumulative masses (as in Figs, la and 2) as a function of time 
for Uranus case U2. Note that the decrease in assumed planetesimai size (and the corresponding reduction in planetesimal random velocities) 
greatly reduces the duration of phase 1. 


the nominal Saturn model (9.8 x 10 6 years) than for the 
nominal Jupiter model (8.0 x 10 6 years). This surprising 
similarity in r end for the nominal Jupiter and Saturn models 
is because the isolation mass is nearly the same for the 
two models (see Section 4). Because the isolation masses 
were chosen to be the same, the value of M z at the conclu- 
sion of the Saturn simulation (17.5 Af©) is close to that for 
the nominal Jupiter case (21.5A/©). 

Case S2 is identical to case SI except that S s is set to 
zero. As in cases J1 and J4, the lifetime of phase 1 and 
the value of Af cross are hardly affected, but the lifetime of 
phase 2 is drastically reduced to 1 x 10 6 years, as compared 
with 7 x 10 6 years for case SI. Thus, the Saturn formation 
time is reduced to only 3.3 X 10° years in this case. 

As illustrated by Figs. 2a and b and discussed above, 
f end depends quite sensitively on cr ^ . Thus, whether Jupiter 
or Saturn first reached runaway gas accretion was deter- 
mined by differences in the actual values of cr imt from the 
ones selected for our nominal models. In principle, we 
might be able to make rough estimates of these differences 
if we knew accurately the high-Z masses of these two 
planets. However, at present there is disagreement as to 
which planet has the larger M z (Zharkov and Gudkova 
1991, Chabrier et al 1992). 

Figure 5a shows the evolutionary history of the “Ura- 
nus” case Ul, where cr init = 0.75 g/cm 2 and a = 19.2 AU. 
Again we simply scaled cr init as a" 2 . The time scale for phase 
1 reaches about 1.5 X 10 7 years, about a factor of 8 longer 
than that for Saturn. Again, M z (r C nd) - 17 A/© and r end = 
2.2 X 10 7 years, a factor of 2.2 longer than that for Saturn. 
Note that Mxy — 1.7 Af© and M z — 12.4M©, comparable 
to the present-day Uranus, after 1.6 x 10 7 years of evolu- 
tion. The period during which Mxy is in the range 1-4 A/© 


lasts 4 million years, from 15 to 19 myr. In a rerun of the 
“Uranus” case with small planetesimals (case U2; Fig. 5b) 
the factor F g reaches a maximum of 7.5 x 10 5 during the 
early phases as compared with a maximum of 1.8 x 10 4 in 
case Ul. As a result, f phl is drastically shortened by a factor 
15, to 1 X 10 6 years. The time for phase 2, however, is not 
much affected, lasting 6.7 x 10 6 and 5.8 x 10 6 years in 
cases Ul and U2, respectively. Also, the values for M cross 
are very similar. However, the time at which Mxy = 1.7 M& 
in case U2 is only 1.6 x 10 6 years, a factor of 10 earlier 
than the comparable time for case Ul. The envelope mass 
stays in the range l-4Af© between 1.5 X 10 6 and 3.5 x 10 6 
years. Thus, a model with characteristics similar to those 
of the present planet can easily be obtained on time 
scales ^ t ncb . The parameters for case U2 give a formation 
time for Uranus that is too short compared with the nomi- 
nal formation times of Jupiter and Saturn; a slightly smaller 
value of crinit would improve the fit. Since the planetesimal 
size has a decisive influence on the time scale for evolution 
of the model for Uranus, that time scale should be consid- 
ered very uncertain, and future work should include con- 
sideration of a range of planetesimal sizes and the evolution 
of the size distribution as a result of accretion. 

4. DISCUSSION 
4.1. Gas Accretion Rate 

Here, we try to understand the factors that control the 
gas accretion rate, especially the conditions that lead to 
runaway gas accretion (phase 3). Despite differences in 
absolute scales, all evolutionary models run to date share 
certain basic characteristics. There are always three phases. 
These phases are distinguished by the temporal behavior 
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of the gas and planetesimal accretion rates, the relative 
magnitudes of these two rates, and the relative magnitudes 
of the cumulative amounts of accreted low- and high-Z 
material. We use these properties, in conjunction with the 
sensitivity of the gas accretion rate to key parameters, to 
infer the factors that control the rate of gas accretion during 
each of these phases. For convenience, we use the terms 
envelope and low-Z mass interchangeably below, as well 
as the terms core and high-Z mass . 

In a formal sense, the rate of gas accretion, as determined 
here, is defined by the new volume of space opened up at 
the outer edge of the planet's envelope by a combination 
of the contraction of constant-mass shells near this edge 
and the expansion of the outer boundary that results from 
the increase in the protoplanet's total mass. The basic 
properties of our evolutionary models can be understood 
in terms of which of these two processes is the dominant 
one and which component of the accretion controls its rate 
of change. During phase 1. the envelope’s mass is small, 
and, except near the end. the planetesimal accretion rate 
is high, always exceeding that of the gas accretion. The 
planet’s mass is increasing rapidly almost solely due to the 
accretion of planetesimals. However, the rate at which 
the outer envelope contracts is greater than the rate of 
expansion of the outer boundary'. Therefore, the rate of 
gas accretion is controlled by the rate at which the envelope 
contracts, which, in turn, is controlled by the energy sup- 
plied by planetesimal accretion. This conclusion is consis- 
tent with the fact that Mxy is more than twice as large at 
the end of phase 1 in case J4 than it is in case Jl. For 
Jupiter models with 100-km planetesimals, Eqs. (14) and 
(15) give f phl * cr- n VF~ g K As the average value of F g in- 
creases weakly with cr imt (because the average mass of the 
planet is larger), the actual calculations are closer to f ph i cc 

& init * ^ 

During phase 2, the envelope s mass is smaller than 
that of the core, the gas accretion rate exceeds that of 
planetesimal accretion, and both rates are nearly constant 
in time and relatively low. This behavior suggests that the 
rate of gas accretion is determined by dR p !dt. The reason 
is that the high-Z material dominates the total mass, while 
the added material is mainly gas. so that the planet’s total 
mass, and therefore R b d , increases only slowly with time. 

To examine the nature of phase 2 more carefully, we 
reran the baseline case but stopped all planetesimal accre- 
tion after certain selected times (1.5, 3.5, and 6.8 myr) 
during phase 2. These runs are denoted as cases Jla, Jib, 
and Jlc, and the masses and accretion rates are plotted as 
a function of time in Figs. 6a and b, respectively. The 
results show that the gas accretion rate (Fig. 6b) jumps 
suddenly by a factor of about 4 in each case right after 
planetesimal accretion is halted. This change is a conse- 
quence of the energy balance within the planet, as the 
planet’s total luminosity remains the same. The length of 


the remaining time in phase 2 is reduced to 1.5, 1.0. and 
0.25 mvr. respectively, in the three cases, corresponding 
in each case to a factor of 4 shorter than the remaining 
time in the baseline case. 

These calculations, along with the previous cases, help 
to define the nature of phase 2. Depending of the principal 
energv source, the length of this phase depends on either 
the gravitational contraction (Kelvin-Helmholtz) time 
scale" of the envelope or the time scale for accretion of 
planetesimals. In the first case, the contraction time is given 
to within a factor of 2 by 

|£ g ravt G.V/ corc .W 

env /iCn 

RL~’ (19) 


where is the gravitational energy of the added mass 
and L is the planet’s luminosity. For the baseline case, in 
the absence of planetesimal accretion. r c ~ 1 mvr. based 
on a radius R = 5 x 10 10 cm. inside of which —909c of the 
envelope mass is contained. This estimate is in agreement 
with the results of case Jla. where the planetesimal accre- 
tion is cut off early in phase 2. The evolution of cases Jib 
and Jlc after cutoff is faster, because less gas remains to 
be accreted to reach iW CTOSS . This time scale is also consistent 
with the phase 2 times of cases J4 (1.3 myr) and S2 (1.0 
myr). in which the energy generation rate from planetesi- 
mals is sharply reduced compared with the standard case. 
The luminosity, which is an important factor in the determi- 
nation of r c . depends on internal properties of the model, 
such as opacity. 

^5 gas is accreted, how'ev^r. the added mass results in 
an increased supply of planetesimals in the feeding zone 
[cf. Eqs. (4) and (7)]. As these accrete onto the protoplanet 
[Eq. (17)], they generate an accretion luminosity, which 
for 5 S = 1 is given by 


L D ,= 


G.V/ core Mz 

Rc ore 


( 20 ) 


The phase 2 luminosities for case Jl, J2, and J3 are in good 
agreement with this expression, given the calculated values 
of M z . However, M z itself is determined by other factors. 
The main physical effect that determines the time scale is 
the mass-luminosity relation that is intrinsic to the struc- 
ture of the protoplanet. The energy loss rate through the 
planet is determined by the rate of radiative transfer, which 
depends on the opacity. In the case of stars with constant 
opacity, standard stellar interiors theory gives L * M 3 
(Clayton 1983, p. 185). In the present case, the situation 
is complicated by the core-envelope structure and by the 
variation of the opacity with temperature and density. Dur- 
ing most of phase 2, where the relevant mass is close to 
the isolation mass, the numerical results are closer to L <*■ 
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FIG. 6. (a) Cumulative masses (as m Fig. la) as a function of time for cases Jla. Jib. and Jlc. The curves extending farthest to the right 

correspond to case Jl. The other curves (left to right) show runs in which planetesimal accretion was arbitrarily stopped at times 1.5. 3.5. and 6.8 
myr, respectively, with all other parameters the same as in case Jl. Crosses show the core mass ( M z ) at these times, (b) Accretion rates as a function 
of time for case Jl and Jla. The curves extending farthest to the right correspond to case Jl. The other curve shows the gas accretion rate in which 
planetesimal accretion was arbitrarily stopped at time 1.5 myr. 


M n , where 4 ^ n ^ 5. In any case, if planetesimal accretion 
supplies most of the radiated luminosity and all high-Z 
material sinks to the core (as in these three cases), the 
time scale is approximately derived as follows: The total 
energy released during phase 2 is 

£«ot - L X r ph2 = G-V/ C ore — , (21) 

‘'core 

where \M Z is the added planetesimal mass during phase 
2. [If contraction of the envelope is an important energy 
source, or dissolved planetesimals do not sink, a relation- 
ship similar to (21) holds provided the masses and radius 
used are adjusted to represent appropriately averaged 
quantities for the planet.] But A M z M\so [cf- Eq. (18)]. 

and Rcorc 00 Mlore f therefore 

f P K2 * (22) 

Contraction of the planet's envelope determines the ther- 
mal energy input for the planet, contributing both directly 
via the release of gravitational potential energy in the enve- 
lope and indirectly because it controls Mxy and thus, via 
Eq. (17), M z and the energy supplied by planetesimal ac- 
cretion. For example, if the luminosity is low, a slow rate 
of contraction is required to bring planetesimals in at a 
sufficient rate to supply this luminosity; therefore, M z also 
is low [Eq. (17)] and the time scale is long. To compare 
with the numerical results, cases J2 and J3 differ in isolation 
mass by a factor of almost 3. Assuming a mass-luminosity 
relation with n = 4§, the duration of phase 2 should go 
approximately as which is reasonably consistent with 


the factor of 43 difference that is actually obtained in these 
two cases. Table V shows a comparison between the above 
estimate for phase 2 time scales (using for definiteness 
the minimum value of luminosity, L min ) and the actual 
numerical results. The similarity of the numerical values 
of the ratio listed in the final column of Table V for runs 
Jl and J5-J7 is a consequence of the energy balance; the 
similarity of the number for runs J2 and J3 supports the 
scaling given by expression (22). Substantially larger ratios 
are obtained for runs J4 and J8 because less energy is 
released by planetesimals if they do not sink to the core, 
and thus, contraction of the envelope supplies most of the 
energy for the planet's luminosity, which implies that Eq. 


TABLE V fl 


run 



imin 


tphl 

■V*/, 3 / (£„.,«,*») 

Jl 

10.0 

11.58 

7.586 x 10~* 

7.S14 x 10 s 

6.97 

3.554 

X 

10-* 

J2 

7.5 

7.52 

6.457 x L0 --# 

4.471 x 10 8 

47.19 

3.004 

X 

1Q-* 

J3 

15.0 

21.28 

1.318 x 10-* 

1.240 x 10 s 

1.06 

3.708 

X 

io- 8 

J4 

10.0 

11.58 

1.413 x 10- r 

4.196 x 10 8 

1.27 

1.048 

X 

L0” 5 

J5 

10.0 

11.58 

6.918 x 10’* 

8.567 x 10 8 

7.02 

3.869 

X 

io- fi 

J6 

10.0 

11.58 

2.512 x 10- 7 

2.360 x 10 8 

2.20 

3.401 

X 

10-* 

J7 

10.0 

11.58 

7.762 x 10~* 

7.636 x 10* 

6.68 

3.624 

X 

10"* 

J8 

10.0 

11.58 

1.238 x 10- 7 

4.601 x 10 8 

0.96 

1.520 

X 

lO’ 5 

SI 

3.0 

11.73 

8.511 x 10-* 

7.115 x 10 8 

7.02 

3.213 

X 

I0" 8 

S2 

3.0 

11.73 

1.862 x 10- 7 

3.252 x 10 8 

1.01 

1.021 

X 

io- 3 

U1 

0.75 

11.92 

8.913 x 10-* 

6.979 x 10 8 

6.70 

3.303 

X 

lO" 6 

U2 

0.75 

11.92 

9.772 x 10“* 

6.365 x 10 8 

5.77 

3.497 

X 

10- 6 


o*iniLZ is in units of g/cm 2 , M is in units of Earth masses, 
L mm is in units of solar luminosity. and r ph2 is in units of myr. 
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(19) rather than expression (22) is the primary determining 
factor for t pb2 . Note that if the contraction rate of the 
planet were to increase suddenly, the rate of accretion of 
solids would also increase, generating an increased plane - 
tesimal accretion luminosity, which would tend to suppress 
the rapid increase in Mxy * Thus, to some extent phase 2 
is self-regulating and stable. 

In all models presented here, the transition to runaway 
gas accretion begins when \f z ~ O.SM^^ anc * Mxy ** 
0.2 A/cross » however, for definiteness and for consistency 
with earlier work, we define the onset of phase 3 at M z = 
Mxy — Across- During the transition and during phase 3 
itself, the rates of gas and planetesimal accretion both 
increase in a quasi-exponential fashion with time, but the 
gas accretion rate grows more rapidly than does the plane- 
tesimal accretion rate. Hence, the cumulative amount of 
low-Z mass in the protoplanet in phase 3 becomes progres- 
sivelv more and more the dominant component [cf. Eqs. 
(16) and (17), but note that in phase 3 when Mxy is ver Y 
large, the accretion rate of pianetesimals is unable to keep 
pace with the expansion of the planet’s accretion zone, so 
Eq. (17) overestimates M z ]. This behavior suggests that 
the gas accretion rate during phase 3 is determined jointly 
by dR P idt and dR bd ldt , with the contraction of R ? being 
the ultimate controlling factor. During late phase 2 and all 
of phase 3, there is an unstable relationship between the 
amount of mass added and the rate at which the envelope 
contracts that leads directly to the quasi-exponential 
growth of the gas accretion rate. Once M z ** Mxy, 
subsequent increase in gas mass produces a comparable 
increase in the protoplanet’s total mass, which, in turn, 
causes a significant increase in R bd - A.s a result, there is 
more gas mass added due to dR?Jdt, which, in turn, leads 
to an Increase in dR p !dt, and so on. During phase 3, the 
PdV work associated with envelope contraction becomes 
the dominant term in the protoplanet’s energy budget, 
whereas the energy associated with planetesimal accretion 
is the dominant term during phase 1 and the largest term 
during phase 2 in most cases considered. The importance 
of the various terms in the energy equation for case J1 is 
shown in Fig. 7. During phase 1 (up to t = 6.3 X 10 5 
years), the radiative loss is almost exactly balanced by 
planetesimal energy deposition, the error is about 2.5%, 
and the gravitational and internal energy terms are negligi- 
ble. During most of phase 2 (up to t ~ 6 x 10 6 years) 
planetesimal deposition stiil represents about two-thirds 
of the energy budget while the typical numerical error 
is 15%. During phase 3, the gravitational energy release 
dominates, and the energy budget cannot be calculated 
accurately because a large amount of new gas is added 
every time step. 

4.2. Comparisons with Other Calculations 

In what ways have our simulations shed new light on 
the formation of the giant planets? What are the key issues 



FIG. ?. Global energy budge: as a function of time for case JL The 
various terms in the energy equauon. in erg/sec. are integrated over the 
entire planetary envelope. Plotted are the radiated luminosity (dashed 
line), the energy deposition from pianetesimals (long- and short-dashed 
line), the rate of change of internal energy {dotted line), the rate of 
change of gravitational energy of the envelope (dot-dashed line), and 
the sum of all contributions (s odd line), which indicates the error. 


that remain to be investigated? In particular, what are the 
critical parameters, the resulting temporal behavior of the 
rate of gas accretion, the time scale to reach runaway gas 
accretion, and the low- and high-Z masses at this point? 

The classical static calculations of Mizuno (1980) and 
the later evolutionary simulations of BP86 were performed 
for solar composition envelopes, grain-dominated opacity 
in the outer envelope, and. of greater importance here, a 
time-invariant rate of planetesimal accretion. The simula- 
tions of this paper are distinguished from those of its prede- 
cessors chiefly by the explicit calculation of the rate of 
planetesimal accretion for situations involving isolated, 
massive protoplanetary embryos surrounded by a swarm of 
much less massive pianetesimals. This choice of accretional 
environment for the for min g giant planets ensures a large 
time variation in the planetesimal accretion rate, as, in 
fact, is well demonstrated by the results of our simulations 
(cf. Figs, la, b). Although our choice, therefore, represents 
an extreme situation, it constitutes one plausible pathway 
by which the giant planets could have formed on time 
scales that are consistent with the lifetime of the solar 
nebula and the age of the Solar System (Lissauer 1987). 

In the older simulations, the gas accretion rate increased 
with time in a quasi-exponential fashion throughout the 
entirety of the simulation, although its rate of increase 
became steeper with increasing time. In this sense, there 
was just one continuous phase of accretion, with the run- 
away portion simply being marked by the time where the 
gas accretion rate first exceeds the planetesimal accretion 
rate. This point lies quite close to that where the low- and 
high-Z masses equal one another. In the current calcula- 
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tions, there are three major phases along the evolutionary 
tracks. During phases 1 and 3. the gas accretion rate in- 
creases in a quasi-exponential fashion with time, but it has 
a nearly constant rate during phase 2. Phase 2 occurs when 
the planet’s feeding zone is almost entirely depleted by 
previous accretion, and its properties depend strongly on 
the mutual interactions of the planetesimal and gas accre- 
tion rates. More generally, allowance for a time-varying 
rate of planetesimal accretion leads to both qualitative and 
quantitative differences between the current calculations 
and their predecessors. 

Runaway gas accretion, however, still begins near the 
point where the low- and high-Z masses are equal. We 
suspect that this similarity arises from a common cause: 
the sensitivity of the gas accretion rate to the location of 
the planet’s outer boundary once the low-Z mass becomes 
a significant fraction of the planet's total mass. Under these 
conditions, there is a strong and unstable relationship be- 
tween the rate of contraction of constant-mass shells in 
the outer envelope and the rate of expansion of the outer 
boundary, leading to runaway gas accretion. 

Previous calculations (BPS6 > showed that the crossover 
mass, A/ cross , was insensitive to the outer boundary parame- 
ters p neb and 7 neb . In the present calculation, M cross is also 
insensitive to p neb (7 neb was not tested). The value that 
BP86 obtained for iV/ cross was modestly sensitive to the 
grain opacity, in the sense that it decreased by 21% when 
the opacity was reduced by a factor of 50. In the present 
case, a factor of 50 reduction in opacity results in practically 
no change in A/ cross , which depends almost solely on M lso 
[Eq. (18)] and. therefore, on a and a [Eq. (14)], but the 
time scale to pass through phase 2 is a factor of 3 shorter. 
In both sets of simulations, the reduction in opacity results 
in an increase in radiated luminosity by a factor of 3, so 
the rate of gravitational contraction must speed up, relative 
to the baseline case, to supply the extra energy. Thus, there 
is a faster rate of gas accretion from the nebula. In our 
simulations, this increase in Xlyy leads to an increase in 
M z [cf. Eq. (17)] so there is no change in A/ cross , but r cross 
decreases significantly, whereas M z is fixed in the simula- 
tions of BP86, so that both and / cross decrease with 

decreasing opacity. The range of values for M CTOSS in BP86 
(ll-29iV/^) is practically the same as the range in the pres- 
ent calculation (10-30 M e ). The values of M aoa computed 
by BP86 also depended on the assumed planetesimal accre- 
tion rate, M z \ a factor of 10 increase in this rate resulted 
in an increase in A/ cross bv 72%. We find, as before, that 
a higher A/ cross is associated with a higher planetesimal 
accretion rate, corresponding in the present case to a 
higher <x inu . 

In the calculations of BP86. the time scale to reach run- 
away gas accretion was determined almost entirely by the 
specified value of M z \ these times ranged from 3 x 10 6 to 
1.12 X 10 8 years. In the current calculations a similarly 


wide range was found, from 2 X 10 6 to 5 X 10 7 years, with 
a strong dependence on cr imt . a, <5 S , and the opacity. This 
time scale is in most cases determined by the length of 
phase 2, which was discussed in Section 4.1. For example, 
for standard Jupiter parameters (cases J1 and J4), a change 
in <5 S from 1 to 0 resulted in a change in time scale from 
8 x 10 6 to 1.6 x 10 6 years. Presumably the “truth” lies 
somewhere inbe tween the extremes of S s = 0 and S s = 1. 
Thus, reasonable parameter choices lead to formation 
times for Jupiter and Saturn of 5 x 10 6 and 7 x 10 6 years, 
respectively. A surprising result of this calculation is the 
similarity of the phase 2 time scales in the standard models 
of Jupiter, Saturn, and Uranus (cases Jl, SI, Ul): they are 
all 7 x 10 6 years. This result is a consequence of the same 
isolation masses being chosen in all three cases, so that 
the model structures at the end of phase 1. the luminosities, 
the gravitational contraction times, and the energy deposi- 
tion by planetesimals in phase 2 were all about the same. 
In retrospect, this is not surprising, as we have already 
shown that the duration of phase 2 is not sensitive to the 
density and temperature of the nebula gas, which are the 
only other variables involved. 

4.3. Implications for the Solar Nebula 

Within the context of our model, it is possible to derive 
tight bounds on the “initial*’ surface density of planetesi- 
mals, cr init , in the outer Solar System. These bounds come 
from the joint constraints of time scale and the current 
high-Z masses of the giant planets. On the one hand. Oi nit 
needs to be smaller than some upper bound or the model 
planets would have accreted more high-Z mass than they 
currently contain. On the other hand, cr ina needs to be 
larger than some lower bound or giant planet formation 
would violate a time scale constraint. This constraint is set 
by the lifetime of the gas component of the solar nebula, 
r sn , estimated from observations to be ^10 7 years. Jupiter 
and Saturn need to reach phase 3 and Uranus and Neptune 
need to reach phase 2 before this time. 

We estimate that <r imt ^ 10 g/cm 2 in the region of the 
solar nebula where Jupiter formed. This value cannot be 
increased by more than about 50% nor decreased by more 
than about 20% without violating one of the two constraints 
cited above. This value is about a factor of 4 larger than 
that given by the minimum-mass solar nebula of Hayashi 
et al. (1985). It is also about a factor of 2 smaller than that 
estimated by Lissauer (1987). based on the assumption 
that Jupiter’s core grew to >1 5A/ 0 before it became iso- 
lated. Our value of cr init is somewhat smaller than that of 
Lissauer because a giant planet can still acquire high-Z 
mass after nearly depleting its feeding zone at the end of 
phase 1 by the outward expansion of its feeding zone during 
phases 2 and 3. As a result, a smaller demand is placed on 
the fully formed Jupiter to scatter the remaining planetesi- 
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mals out of the Solar System. Put another way, Jupiter 
should have moved only a small fraction of its initial dis- 
tance from the Sun due to this clearing out of planetesimals, 
as its total mass would have been significantly larger than 
the cumulative mass of removed planetesimals. 

Since the terrestrial planets have relatively low masses 
and are more deeply embedded in the Sun*s gravitational 
well than is Jupiter, they would not have been able to 
effectively scatter planetesimals out of the region of the 
inner Solar System. Thus, cr in the inner Solar System 
should have been dose to the values provided by the mini- 
mum-mass solar nebula (e.g.. Lissauer 1987). Therefore, 
cr varied only slowly with distance from the Sun throughout 
the inner Solar System and out to the distance of Jupiter, 
aside from a possible discontinuity caused by the condensa- 
tion of water ice near the vicinity of Jupiter (Lissauer 1987). 
However, our two basic constraints on cr in the region 
where the giant planets formed imply that cr oc a' 2 in the 
outer solar nebula if the high-Z mass is actually the same 
for Jupiter and Saturn. Thus, the values of cr given by our 
simulations tend to approach those of the minimum-mass 
solar nebula of Hayashi et cil (1985) with increasing dis- 
tance from the Sun; i.e., the enhancement factor was largest 
for Jupiter. 

Our inferred dependence of cr on distance from the Sun 
in the giant planet region ( a~ z ) is much steeper than that 
expected for a fully viscously evolved disk (no steeper than 
a' 0 - 5 ; Ruden and Lin 1986) and is more nearly comparable 
to the distribution in a disk that has just formed by collapse 
from a rotating molecular cloud core (cT //4 ; Cassen and 
Moosman 1981) or that inferred from the spectral energy 
distributions of young stars (a ~' 2 ; Beckwith et ai 1990). 
This comparison suggests either that little viscous evolution 
occurred in the outer region of the solar nebula or that 
the radial evolution of solids became quickly decoupled 
from that of gases in this region of the solar nebula. By 
contrast, the relatively flat profile of cr inferred for the 
inner region of the solar nebula implies just the opposite. 
The viscous evolution time scale can increase with distance 
from the Sun (Ruden and Lin 1986), so such a surface 
density profile may be realistic. We add the strong caveat 
that the conclusions drawn here and elsewhere in this paper 
depend strongly on the validity of the basic assumptions of 
our simulations: isolated embryos and no systematic radial 
motions of planetesimals. 

Finally, we examine more carefully the conditions under 
which our scenario for giant planet formation satisfies cur- 
rent estimates of f sn . Our nominal models of Jupiter and 
Saturn reach runaway gas accretion on time scales of 8.0 X 
10 6 and 9.8 X 10 6 years, respectively. We set S s = 1 in these 
simulations; i.e., we assume that vaporized planetesimal 
material sank to the core interface. The opposite extreme, 
that of no sinking, shortens the time to reach runaway gas 
accretion by about a factor of 4 (cf. Fig. 3). This time 


interval can also be shortened by choosing a slightly larger 
value of <r imt . Thus, it is relatively easy to find models for 
Jupiter and Saturn that reach phase 3 in a time interval 
comparable to or less than r ra . Note, however, that these 
time scales again depend on our basic assumptions. For 
example, if the planetesimal random velocities were as- 
sumed to be large, the time scale for phase 1 could well 
exceed that for phase 2. 

At a time of 1.6 x 10 7 years, early in phase 2, our nominal 
model reaches a point where its high- and low-Z masses 
are comparable to those of present Uranus. With smaller 
planetesimals this time w'as reduced to 1.5 x 10 6 years! 
Note that S % does not have much influence on the time 
scale in this case. Thus, our formation scenario for Uranus 
is compatible with current estimates of f sn . Neptune has 
always been the most difficult planet to form quickly 
enough because of the steep increase in the formation time 
scale with increasing distance from the Sun. Comparison 
of our nominal models of Jupiter, Saturn, and Uranus (cf. 
Figs. la. 4, 5a) indicates that the longevity of phase 1. r ph i, 
scales roughly as of. Since Neptune has high- and low- 
Z masses comparable to those of Uranus (Zharkov and 
Gudkova 1991). we expect that this scaling law will give 
us an approximate formation time for Neptune. Extrapola- 
tion of our results suggests that f P hi 555 3.7 X 10' years for 
Neptune for large (100-km) planetesimals but only 3.7 x 
10 6 years for small (1-km i planetesimals. If accretion was 
dominated by reasonably small planetesimals. we get a 
time scale for Neptune that is compatible with t sn . Because 
phase 2 lasts for a long period, these calculations partially 
resolve the issue of why Lranus and Neptune ha^e very 
similar properties. However, the observed similarity be- 
tween Mz values of Jupiter. Saturn, Uranus, and Neptune 
must be considered a consequence of the ^initial” distribu- 
tion of condensed mass within the protoplanetary disk 
rather than a result of fundamental growth processes of 
the giant planets. 

4.4. Implications for Giant Planets 

As illustrated in Fig. le. the capture radius for planetesi- 
mals starts to significantly exceed the core radius when the 
accreted mass of high-Z material reaches ** 2 Shortly 
afterward, at a core mass M<ns *** 2.6Af$, planetesimals 
dissolve in the envelope. The precise value of M& s de- 
creases modestly as the planetesimal size decreases, and 
it is smaller by a few tens of percent in case SI as compared 
with case Jl. Overall, its range is 2-4 A/a. This finding 
strongly suggests that much of the high-Z mass of the 
giant planets is located (although not necessarily uniformly 
distributed) within their envelopes, rather than in a segre- 
gated core. This prediction is in good accord with the 
results of recent interior models (Zharkov and Gudkova 
1991, Chabrier et ai 1992). which yield interior core masses 
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FIG. 8. Rate of energy deposition by planetesiraals in the planetary envelope as a function of the distance from the planet’s center. R. at a 
time of 3.5 x 10 5 years in the evolution of case Jl. The planet’s radius R p — 4 x 10 l ° cm. Energies are given in ergs/g/sec. Plotted are the total 
energy (solid line), the latent heats (dotted lines) delivered by organics (ecHos) and absorbed by H : 0 and rock, the contribution from gas drag 
(dashed line), the gravitational potential energy (long- and short-dashed line), and the kinetic energy (dot-dashed tine). The core mass is 1.4 
and the envelope mass is 8.6 x 10 _o \U. (a) The curves are overplotted using the same scale, (b) The same curves are displaced vertically by 
arbitrary amounts for clanty. The number at the upper end of each curve gives the logarithm of the maximum rate of energy deposition in units 
of erg/g/sec. 


in the range 1 to 8 Me for Jupiter and Saturn. The core 
mass could be augmented over M dis by some settling of 
material to the core or by the occasional accretion of a 
very massive planetesimal; nevertheless, the “truth” prob- 
ably lies closer to S $ = 0 than S s = 1. This point is illustrated 
in Figs. 8 and 9, which show where in the protoplanet most 
of the accretion energy of the planetesimals is deposited. 
At an early time, during phase 1 when M z = IM 9 , most 
of the energy is deposited at the core boundary (Fig. 8 ), 
while at a later time during phase 2 much of the energy is 
deposited at 10 core radii (Fig. 9). 

It is likely that some of the dissolved high-Z material 
ultimately gets mixed throughout the planets’ envelopes. 
The models of the present study have extensive convection 
zones, similar to those reported by BP 86 . Furthermore, 
the calculations of the subsequent evolution (BP 86 ) show 
that the envelopes become fully convective during their 
contraction period after accretion has ceased. Thus, an- 
other consequence of our simulations is that elements that 
were derived primarily from the accretion of planetesimals 
(e.g., carbon) should have abundances in the atmospheres 
of the giant planets with respect to H that exceed solar. 
There is good evidence that this is the case (Gautier and 
Owen 1989). Furthermore, the variation in the degree of 
enhancement above solar proportions among the four giant 
planets is qualitatively and semiquantitatively in accord 
with the expectations of planetesimal dissolution (Pollack 


et aL 1986, Podolak et ai 1988. Simonelli et al 1989). Note, 
however, that the present calculations do not account for 
the redistribution of heavy elements through the envelope, 
nor are they carried to the point where the mass approaches 
that of Jupiter; therefore, a detailed comparison of these 
results with observations of heavy element abundances in 
the giant planets is not warranted. Regarding the mixing, 
a possibly significant factor is the radial gradients in mean 
molecular weight induced by deposition of heavy elements, 
an effect that may limit the luminosities of present-day 
Uranus and Neptune (Hubbard et al 1995). We intend to 
account for mixing of the heavy elements, including the 
effect of composition gradients, using more self-consistent 
calculations in the future. 

Zharkov and Gudkova (1991) suggested that the rela- 
tively low core masses of Jupiter and Saturn derived from 
their interior modeling were inconsistent with the critical 
high-Z masses of 10-15 A /9 predicted by core instability 
formation models. Therefore, they concluded that Jupiter 
and Saturn did not obtain their low-Z masses by means of 
the runaway gas accretion process that occurred once the 
critical core mass was attained. In the present calculation, 
all high-Z mass was assumed to end up in the core, although 
the energy deposition in the envelope was taken into ac- 
count. Future calculations including mass deposition in the 
envelope are necessary to test Zharkov and Gudkova’s 
conclusions. 
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FIG. 9. Same as Fie. 8 except that r = 4.0 x 10* years. * p = 5.076 x 10". the core mass is 13.2 W,. and the envelope mass is ^V/^ Note that 
in this case a larger fraction of planetesimal energy is deposited well above the core because the planetesimais do not reach the core inta . 


5. CONCLUSIONS 

In the evolution of a giant planet, allowance for a vari- 
able rate of planetesimal accretion (M z ) results in an im- 
portant qualitative difference from the case in which M z 
is assumed to be constant. In the case of constant M z , 
there are two phases. In the first phase, M z dominates the 
energy generation and M z > Mxy- In second phase, 
M z < Mxy . rapid gas accretion occurs, and gravitational 
contraction dominates the energy production. In the new 
scenario, there are three phases. In the first, M z dominates 
the energy production and increases rapidly to a maximum, 
then declines as the isolation mass is reached. In the sec- 
ond, M z > Mxy , M z declines to a low and nearly constant 
level, Mxy ** 2-4 M z , but M z still dominates the energy 
production. The third phase is analogous to the second 
phase of the previous scenario (BP86) with constant M z 
and sets in, by definition, when M z = Mxy = Across - Most 
of the evolution is generally spent in phase 2. Our results 
suggest a plausible “pathway" by which Jupiter and Saturn 
could reach phase 3 and undergo rapid gas accretion after 
a total elapsed time of a few million years, shorter than 
the lifetime of a solar nebula. We can explain the inferred 
bulk composition of Uranus and Neptune as the result of 
the dissipation of the gas component of the solar nebula 
while the planets were still in the long-lived phase 2, with 
still relatively small Mxy • Our simulations have also con- 
firmed that a key condition for the inception of runaway 
gas accretion in the parameter regime explored is that the 
mass of the low-Z envelope is a substantial fraction of 


the high-Z core. The feedback process responsible for the 
runaway is the coupling of the contraction of the material 
of the envelope and the expansion of the outer boundary 
of the envelope (the accretion radius) with Mxy- Finally, 
the results of our simulations, which take into account the 
dissolution of planetesimals in the gaseous envelope, are 
consistent with the partitioning of high-Z mass between 
the core and the envelope in the current giant planets 
and with the supersolar abundance ratios of planetesimal- 
derived elements in the atmospheres. 

The actual rate at which the giant planets accreted small 
planetesimals is probably intermediate between the con- 
stant rates assumed in most previous studies and the highly 
variable rates found in the present study. The main assump- 
tions in the accretion model are (1) an isolated embryo, 
(2) an initial phase of runaway accretion of solids, (3) small 
random velocities and sizes of planetesimals, and (4) no 
planetesimal migration into or out of the current feeding 
zone. Given this model, the simulations provide strong 
and interesting constraints on the initial surface density of 
planetesimals in the outer solar nebula; for our standard 
parameters we find that oj apiter ~ 10 g cm ", asatum ~ ^ g 
cm' 2 , and <r Uranils ~ 0.75 g cm' 2 . The corresponding forma- 
tion times (assuming a planetesimal radius — 100 km) are 
8 x 10 6 , 1 X 10 7 , and 1.6 x 10 7 years. The high-Z masses 
and the evolutionary times are extremely sensitive to cr mit . 
In our opinion, these constraints should not be taken too 
literally, given the specific set of assumptions under which 
the calculations were performed. The above formation 
times should be regarded as conservatively long. For exam- 
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pie, if our assumption that all accreted high-Z material 
eventually falls to the core is relaxed, and material is al- 
lowed to remain in the envelope (<5 S = 0), formation times 
for the giant planets can be reduced to 3 X 10° and 3.5 x 
10 6 years for Jupiter and Saturn, respectively. In the case 
of Uranus, the formation time is dominated by the length 
of phase 1 rather than phase 2, and S s has little effect. If, 
however, the assumed planetesimal size is reduced by a 
factor of 100, the time for the formation of Uranus is 
reduced by a factor of 10. but the times for Jupiter and 
Saturn are little affected. Or or ^ could be increased 
slightly, resulting in reduced formation times for all planets, 
but it cannot be increased by more than 50%; otherwise, 
the high-Z masses become unacceptably large. On the 
other hand, allowing for planetesimal migration will slow 
growth rates for fixed M z total, as more planetesimals are 
accreted at later epochs. There is a strong need to perform 
follow-up calculations in which 1 1 ) alternative planetesimal 
accretion scenarios are examined, such as the allowance 
for stirring by nearby competing embryos, (2) the impact 
of dissolved planetesimal material on the properties of the 
envelope, such as the opacity and convective instability, 
are investigated, (3) the effect of occasional impacts of very 
massive planetesimals are assessed, and (4) hydrodynamic 
effects in the envelope evolution are examined. 
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